{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Laminar pipe flow - Hagen–Poiseuille solution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this tutorial we use icoFoam to compute a laminar pipe flow, then we compare the numerical solution with the analytical solution. <br>\n",
    "\n",
    "You will find the instructions of how to run this case in the file README.FIRST. <br>\n",
    "\n",
    "By the way, to plot the numerical solution you need to use the utility sample to get the solution at the end of the pipe."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Remember, to execute a cell in ipython you need to press shift-enter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The %matplotlib magic is used to specify which matplotlib backend we want to use. \n",
    "Most of the time, in the Notebook, you will want to use the inline backend which will embed plots inside the Notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "ModuleNotFoundError",
     "evalue": "No module named 'matplotlib'",
     "output_type": "error",
     "traceback": [
      "\u001b[31m---------------------------------------------------------------------------\u001b[39m",
      "\u001b[31mModuleNotFoundError\u001b[39m                       Traceback (most recent call last)",
      "\u001b[36mCell\u001b[39m\u001b[36m \u001b[39m\u001b[32mIn[1]\u001b[39m\u001b[32m, line 1\u001b[39m\n\u001b[32m----> \u001b[39m\u001b[32m1\u001b[39m \u001b[43mget_ipython\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[43m.\u001b[49m\u001b[43mrun_line_magic\u001b[49m\u001b[43m(\u001b[49m\u001b[33;43m'\u001b[39;49m\u001b[33;43mmatplotlib\u001b[39;49m\u001b[33;43m'\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[33;43m'\u001b[39;49m\u001b[33;43minline\u001b[39;49m\u001b[33;43m'\u001b[39;49m\u001b[43m)\u001b[49m\n",
      "\u001b[36mFile \u001b[39m\u001b[32m~/miniconda3/envs/tttt_playground/lib/python3.12/site-packages/IPython/core/interactiveshell.py:2504\u001b[39m, in \u001b[36mInteractiveShell.run_line_magic\u001b[39m\u001b[34m(self, magic_name, line, _stack_depth)\u001b[39m\n\u001b[32m   2502\u001b[39m     kwargs[\u001b[33m'\u001b[39m\u001b[33mlocal_ns\u001b[39m\u001b[33m'\u001b[39m] = \u001b[38;5;28mself\u001b[39m.get_local_scope(stack_depth)\n\u001b[32m   2503\u001b[39m \u001b[38;5;28;01mwith\u001b[39;00m \u001b[38;5;28mself\u001b[39m.builtin_trap:\n\u001b[32m-> \u001b[39m\u001b[32m2504\u001b[39m     result = \u001b[43mfn\u001b[49m\u001b[43m(\u001b[49m\u001b[43m*\u001b[49m\u001b[43margs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43m*\u001b[49m\u001b[43m*\u001b[49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n\u001b[32m   2506\u001b[39m \u001b[38;5;66;03m# The code below prevents the output from being displayed\u001b[39;00m\n\u001b[32m   2507\u001b[39m \u001b[38;5;66;03m# when using magics with decorator @output_can_be_silenced\u001b[39;00m\n\u001b[32m   2508\u001b[39m \u001b[38;5;66;03m# when the last Python token in the expression is a ';'.\u001b[39;00m\n\u001b[32m   2509\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mgetattr\u001b[39m(fn, magic.MAGIC_OUTPUT_CAN_BE_SILENCED, \u001b[38;5;28;01mFalse\u001b[39;00m):\n",
      "\u001b[36mFile \u001b[39m\u001b[32m~/miniconda3/envs/tttt_playground/lib/python3.12/site-packages/IPython/core/magics/pylab.py:103\u001b[39m, in \u001b[36mPylabMagics.matplotlib\u001b[39m\u001b[34m(self, line)\u001b[39m\n\u001b[32m     98\u001b[39m     \u001b[38;5;28mprint\u001b[39m(\n\u001b[32m     99\u001b[39m         \u001b[33m\"\u001b[39m\u001b[33mAvailable matplotlib backends: \u001b[39m\u001b[38;5;132;01m%s\u001b[39;00m\u001b[33m\"\u001b[39m\n\u001b[32m    100\u001b[39m         % _list_matplotlib_backends_and_gui_loops()\n\u001b[32m    101\u001b[39m     )\n\u001b[32m    102\u001b[39m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[32m--> \u001b[39m\u001b[32m103\u001b[39m     gui, backend = \u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43mshell\u001b[49m\u001b[43m.\u001b[49m\u001b[43menable_matplotlib\u001b[49m\u001b[43m(\u001b[49m\u001b[43margs\u001b[49m\u001b[43m.\u001b[49m\u001b[43mgui\u001b[49m\u001b[43m)\u001b[49m\n\u001b[32m    104\u001b[39m     \u001b[38;5;28mself\u001b[39m._show_matplotlib_backend(args.gui, backend)\n",
      "\u001b[36mFile \u001b[39m\u001b[32m~/miniconda3/envs/tttt_playground/lib/python3.12/site-packages/IPython/core/interactiveshell.py:3787\u001b[39m, in \u001b[36mInteractiveShell.enable_matplotlib\u001b[39m\u001b[34m(self, gui)\u001b[39m\n\u001b[32m   3784\u001b[39m     \u001b[38;5;28;01mimport\u001b[39;00m \u001b[34;01mmatplotlib_inline\u001b[39;00m\u001b[34;01m.\u001b[39;00m\u001b[34;01mbackend_inline\u001b[39;00m\n\u001b[32m   3786\u001b[39m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[34;01mIPython\u001b[39;00m\u001b[34;01m.\u001b[39;00m\u001b[34;01mcore\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m pylabtools \u001b[38;5;28;01mas\u001b[39;00m pt\n\u001b[32m-> \u001b[39m\u001b[32m3787\u001b[39m gui, backend = \u001b[43mpt\u001b[49m\u001b[43m.\u001b[49m\u001b[43mfind_gui_and_backend\u001b[49m\u001b[43m(\u001b[49m\u001b[43mgui\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43mpylab_gui_select\u001b[49m\u001b[43m)\u001b[49m\n\u001b[32m   3789\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m gui != \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[32m   3790\u001b[39m     \u001b[38;5;66;03m# If we have our first gui selection, store it\u001b[39;00m\n\u001b[32m   3791\u001b[39m     \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mself\u001b[39m.pylab_gui_select \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n",
      "\u001b[36mFile \u001b[39m\u001b[32m~/miniconda3/envs/tttt_playground/lib/python3.12/site-packages/IPython/core/pylabtools.py:338\u001b[39m, in \u001b[36mfind_gui_and_backend\u001b[39m\u001b[34m(gui, gui_select)\u001b[39m\n\u001b[32m    321\u001b[39m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[34mfind_gui_and_backend\u001b[39m(gui=\u001b[38;5;28;01mNone\u001b[39;00m, gui_select=\u001b[38;5;28;01mNone\u001b[39;00m):\n\u001b[32m    322\u001b[39m \u001b[38;5;250m    \u001b[39m\u001b[33;03m\"\"\"Given a gui string return the gui and mpl backend.\u001b[39;00m\n\u001b[32m    323\u001b[39m \n\u001b[32m    324\u001b[39m \u001b[33;03m    Parameters\u001b[39;00m\n\u001b[32m   (...)\u001b[39m\u001b[32m    335\u001b[39m \u001b[33;03m    'WXAgg','Qt4Agg','module://matplotlib_inline.backend_inline','agg').\u001b[39;00m\n\u001b[32m    336\u001b[39m \u001b[33;03m    \"\"\"\u001b[39;00m\n\u001b[32m--> \u001b[39m\u001b[32m338\u001b[39m     \u001b[38;5;28;01mimport\u001b[39;00m \u001b[34;01mmatplotlib\u001b[39;00m\n\u001b[32m    340\u001b[39m     \u001b[38;5;28;01mif\u001b[39;00m _matplotlib_manages_backends():\n\u001b[32m    341\u001b[39m         backend_registry = matplotlib.backends.registry.backend_registry\n",
      "\u001b[31mModuleNotFoundError\u001b[39m: No module named 'matplotlib'"
     ]
    }
   ],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now import the libraries numpy and matplotlib. <br>\n",
    "NumPy is the fundamental package for scientific computing. <br>\n",
    "Matplotlib is used for plotting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now use loadtxt to read the input file.  <br>\n",
    "We type np.loadtxt because loadtxt belongs to numpy which we imported with the name np.<br>\n",
    "The skiprows option is used to skip header lines.  In this case, the input file does not contain a header.<br>\n",
    "All the indormation is saved in the variable data, which is a numpy array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "data = np.loadtxt('../postProcessing/sampleDict/20/s2_U.xy', skiprows=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now compute the velocity profile according to Hagen–Poiseuille solution. <br>\n",
    "\n",
    "$$v_{radial} = v_{max} \\left[ 1 - \\dfrac{r^2}{R_{max}^2} \\right] $$\n",
    "\n",
    "where $v_{max}$ is the maximum axial velocity obtained from the simulation, $R_{max}$ is the maximum radius, and $r$ is the local radius."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "vmax = 1.4\n",
    "rmax = 0.5\n",
    "\n",
    "x = np.linspace(-1*rmax, rmax, 100)\n",
    "sol = vmax*(1 - x**2/rmax**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now plot a sample velocity profile.  <br>\n",
    "To do so, we use the plot function. <br>\n",
    "As this function belongs to matplotlib.pyplot, we append plt to plot as follows: plt.plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fef373c8110>]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEACAYAAABWLgY0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XucleP+//HXp5MichZREhUpx2InNeSQYzknskVyzHb4Ufb3S7MP369tsym0S0oph6jQAQm1kB1ld5AOSnaJTeIrpKKm6/fHtWL2mFprZu61rnWv9X4+HvMwq7lb876t5jPX+tzXdd3mnENERPJTtdABREQkc1TkRUTymIq8iEgeU5EXEcljKvIiInlMRV5EJI+lLPJmNszMVpnZ+ymOa21mG83s3OjiiYhIVaQzkh8OnLqtA8ysGvAX4JUoQomISDRSFnnn3HTgmxSH9QbGAl9GEUpERKJR5Z68me0DdHHODQKs6pFERCQqUVx47Q/0KfVYhV5EJEfUiOA5jgZGm5kBuwOnmdlG59yEsgeamTbKERGpBOdcpQbQ6Y7kja2M0J1zByQ/GuP78teVV+BLHZ+3H/369QueoVDPb8MGx/jxjt/+1rHrro7DDnPccotj4kTH6tXRnN/33zsSCUe/fo7jj3fsuKOjc2fH4487vvkm/P+DuL52Or/UH1WRciRvZk8BRcBuZvYJ0A+o5eu1G1K2hlcpjUgFzZ8Pw4bBk0/CwQfD+efDn/4E++0X/feqWxc6dPAfxcXwzTcwaRKMGwc33ghnnglXXum/Xk0rUCRHpCzyzrlu6T6Zc+6KqsURSW3zZnjpJbjvPvjoI+jRA959Fw44ILs5dtkFunf3H199BU88Ab17+3y33gqXXAK1a2c3k0hZGm9EqKioKHSEjAp9fps3w+jRcOihcNdd0KsX/OtffuQeRYGvyvntvjvcdJN/Z/HQQzB2LDRuDPffD+vXVz1bVYV+7TIt38+vKqyq/Z4KfTMzl83vJ/nBOZg4Ef77v6FOHV/UTz4ZLMfncb3/PvTrBzNnwn/9F1x1FdSsGTqVxJGZ4Sp54VVFXnLaggXwu9/B55/DX/7i+965XtzLmjUL7rjDn0P//v4XlEhFqMhL3lm7Fu68019QvfNOuPZaqBHFhN9AnIPx432vvlUrePhhaNAgdCqJi6oUefXkJedMmQItW/rZKwsX+ouZcS7w4N99dOni35kcdhgcfjgMGeKvM4hkkkbykjPWrvUXL197DR55BE7d5rZ48fbBB366Zd26MGJEZqZ8Sv7QSF5ib9YsOOIIP7KdPz+/Czz4GUL/+Ad07AhHHw1jxoROJPlKI3kJyjk/zfCvf4WBA/1ipkIza5afU9+hg59+qbn1UpYuvEosffcdXHEFfPKJn1fesGHoROF8/z307OkXd22ZYy+yhdo1EjuLF0Pr1rDHHvDWW4Vd4AF23NEv9OreHY49Fl59NXQiyRcq8pJ1r70G7dtDnz4waBBst13oRLnBzF94HjPGF/u//z10IskHatdIVg0e7Df3euYZ34OW8i1bBmed5S/MPvBA/KeQStWoJy85zzm/tH/cOL+5WJMmoRPlvm+/hQsugO23h6ef9ls6SGFST15y2qZN/qLi66/D22+rwKerXj2/lfEOO8App/jFYSIVpSIvGbVhA5x3Hvz73zB1qt+tUdJXqxaMGuXn0rdvD198ETqRxI2KvGTMunXQubO/sDp+vB+RSsVVq+bXElx4ob+O8dlnoRNJnOhyjmTE2rVw9tmwzz5+2b4uHFaNmd+orXZtP6KfOhUaNQqdSuJAP3oSuR9+gDPO8L33Rx+F6tVDJ8oft93mC31REbzxhtYXSGoq8hKp9et9i6ZxYxg6VPc6zYTevaGkxE+vfOMN/25JZGtU5CUyP/3k957ZYw9/c20V+My56SZ/UXtLod9zz9CJJFepyEskSkrg0kv9bJCRI9WiyYa+ff07p1NO8YW+Xr3QiSQXaTGUVJlzcMMN/gYfL7+sXRSzyTnfvvngA5g8Wf/v85VWvEpQf/wjPP+8H03utFPoNIWnpAS6dfOLzp59Vu+i8lFGV7ya2TAzW2Vm72/l693MbF7yY7qZtaxMEImnYcPg8cf9CF4FPozq1X2L7Ntv4cYb/eheZIt0Lo0NB7Z1n56PgfbOucOAPwOPRhFMct9rr8Hvf+/3oqlfP3Sawrbddn5foDfegAEDQqeRXJLywqtzbrqZbXXZhXPunVIP3wF0D/oCsHChbxGMGQPNmoVOI+AvvL74IrRt66ewdu4cOpHkgqgnufUEXo74OSXHrF4NZ54Jf/ubtgvONY0awQsv+A3hZs8OnUZyQWRTKM3sBKAH0G5bxxUXF//8eVFREUVFRVFFkCzYMhe+a1d/YwvJPa1b+5uxnHMOzJwJe+0VOpFUVCKRIJFIRPJcac2uSbZrJjrnWm3l662AcUAn59yybTyPZtfE3HXXwcqVfsMxLXbKbXfd5fe4mTrVr1+Q+MrGfvKW/CjvmzfEF/ju2yrwEn+PPAKJBDz5pAp8HBQX+62db7ghdBIJKeVI3syeAoqA3YBVQD+gFuCcc0PM7FHgXGAF/hfBRudcm608l0byMfXOO35XyenToWnT0GkkXd9/D8ccA7feCldeGTqNVJYWQ0lGrV4NRx0FDz2kGRtxtHgxHH88vPIKHHlk6DRSGbr9n2TMltWUl1yiAh9XzZvDwIH+gvn//V/oNJJtGsnLNt11l2/RTJmiG3/E3c03w9KlMGGCrqnEjdo1khGvvw6XXebnW2saXvxt3OjvKnX++b5HL/GhIi+RW7XK928ffxxOOil0GonK8uXQpg1MmuT/K/GgnrxEavNmP4K//HIV+Hyz//5+oVTXrn5DM8l/GsnLr9x7r1/slEioD5+vrrvOX4R9+ml/k3DJbWrXSGTmzPF3Gpo1y4/6JD+tXw9HHw133OHv6CW5TUVeIrF+vZ8P//vf6we/EMydCyefrF/ocaAiL5G48Ub48ku9hS8k994LEyfCtGm6o1Qu04VXqbIpU/wWtYMGqcAXkltu8cX9vvtCJ5FM0Uhe+PZbaNnS38rv5JNDp5FsW7HC9+cTCWjRInQaKY/aNVIlV17pt6IdNCh0Egnl0Uf9LqMzZkDNmqHTSFlq10ilvfSS32/8r38NnURC6tkTdtsN7rkndBKJmkbyBWzNGt+mGTkSTjghdBoJbeVKP7vqtdegVbm3B5JQ1K6RSrnqKt+mGTgwdBLJFcOGweDBvm2jhXC5Q+0aqbBp0/z+4nffHTqJ5JIrroCddoIBA0InkahoJF+A1q3zb8f794czzwydRnLNsmX+blLvvgtNmoROI6B2jVTQ7bf7/uvTT4dOIrnq3nth8mTfn9e6ifBU5CVtc+fCqafC/Pmw556h00iu2rTJj+ZvvBF++9vQaURFXtJSUgJt20KvXrqps6T23nu+nbdggZ9eKeHowqukZcgQP5umR4/QSSQOjj4aLrwQ+vYNnUSqQiP5AvHFF35O/LRpcOihodNIXHz7LRxyCDzzDLRrFzpN4dJIXlK69VbfolGBl4qoVw8eeACuucbfI1biJ2WRN7NhZrbKzN7fxjEPmtlSM5trZodHG1GqKpGA6dPhzjtDJ5E4uuAC2GcfePjh0EmkMtIZyQ8HTt3aF83sNKCJc+4g4GpgcETZJAIbN0Lv3nD//bDDDqHTSByZwUMPwf/8j2/7SbykLPLOuenAN9s4pDMwMnnsu0A9M9srmnhSVX//O9SvD+eeGzqJxFmzZr7dd/vtoZNIRUWxO0UDYGWpx58l/2xVBM8tVbBqFfz5z/Dmm1rQIlV3553QvLlv/ekibHxkfQui4uLinz8vKiqiqKgo2xEKxh13wOWXw8EHh04i+aBuXX8HqRtugH/+U7cLzKREIkEikYjkudKaQmlmjYCJzrlfbUBqZoOBac65Z5KPFwMdnHO/GslrCmX2vPcenH02LF7sN5wSiYJz0L69XwXbs2foNIUjG1MoLflRngnAZckgxwJryivwkj3OwU03wZ/+pAIv0TLzG9vdeSd8913oNJKOdKZQPgX8A2hqZp+YWQ8zu9rMegE4514C/mVmHwGPANdlNLGk9Oyz8MMPvlUjErWjjoLTTvOzbST3acVrnlm/3l8cGzXKv60WyYTPP/crqLUdcXZoxav87P77oXVrFXjJrL339quob7stdBJJRSP5PLJqld9nZOZMja4k87a8a3zySU2pzDRtNSwAXHcdbLed32tEJBtGjfL3CJ4xQ2sxMklFXli8GI4/3v9Xe39Ltmze7Lck7tvXb0ssmaEiL3Tu7N8yq0cq2TZ1Klx1FSxc6N9JSvR04bXAvfUWzJvnNyITybYTT/SrqgcNCp1EyqORfMw5B8cdB9deC927h04jhWr+fDjpJFi6VAvwMkEj+QI2YQKsXQvduoVOIoWsZUt/g/i//S10EilLI/kYKymBVq3gnnv8DZdFQlq+3K+GXbgQ9tJm45HSSL5AjRoFu+4KZ5wROokI7L8/XHqptjvINRrJx9SPP0LTpvDUU74nL5ILvvzSX4R97z1o3Dh0mvyhkXwBGjLE90FV4CWX7LknXH89/PGPoZPIFhrJx9C6dXDggfDii3DEEaHTiPynNWvgoIP8HaSaNQudJj9oJF9gBg70I3gVeMlFO+8Mt9wCpW4CJwFpJB8z333nR/GJhN+MTCQXrV3r/51OmeJngEnVaCRfQAYM8PORVeAll9WtC336QL9+oZOIRvIxsmaNHx3NmOF7niK5bP16/+91wgQ/f14qTyP5AjFggF/0pAIvcVCnjt+dUjNtwtJIPia2zFiYMcOPjkTiYMMGfwObSZM0UaAqNJIvAA8+6Fe2qsBLnNSu7Xvzf/hD6CSFSyP5GFAvXuJsS29eo/nK00g+z20ZxavASxzVqaPRfEgayee477+HAw6At9/2e9WIxNH69b43P3my5s1XRsZH8mbWycwWm9kSM+tTztd3MrMJZjbXzOab2eWVCSO/NmiQvxmDCrzEWZ06fhXs//5v6CSFJ+VI3syqAUuAjsC/gVlAV+fc4lLH3AHs5Jy7w8x2Bz4E9nLObSrzXBrJV8C6dX4U/+qrfjMykThbu9b/e37rLe1pU1GZHsm3AZY651Y45zYCo4HOZY5xwI7Jz3cEvi5b4KXihg6F3/xGBV7yQ926/j7Ed98dOklhqZHGMQ2AlaUef4ov/KU9DEwws38DdYGLoolXuH78Ee69F55/PnQSkej07u1788uX+5uMSOalU+TTcSowxzl3opk1AV41s1bOubVlDywutTVdUVERRUVFEUXILyNHwqGHwtFHh04iEp2dd4ZrrvG3rBw0KHSa3JVIJEgkEpE8Vzo9+WOBYudcp+TjvoBzzt1T6phJwN3OubeTj18H+jjn3ivzXOrJp6GkBJo3h2HDoH370GlEorV6te/JL1wI9euHThMPme7JzwIONLNGZlYL6ApMKHPMCuCkZJi9gKbAx5UJJPDcc/4OO8cfHzqJSPT22AMuuQT69w+dpDCkNU/ezDoBA/C/FIY55/5iZlfjR/RDzGxvYASwd/Kv3O2ce7qc59FIPgXnfIumuBjOOit0GpHMWL7c70z58cdQr17oNLmvKiN5LYbKMa++CjffDO+/D9W0HlnyWPfu0KKF36lStk1FPo907AiXX+5/AETy2fz5cMopfjRfp07oNLlNe9fkiVmz4KOPoGvX0ElEMq9lS9+aHDEidJL8ppF8DrnwQmjbFm66KXQSkex46y244gpYvBiqVw+dJndpJJ8HPv4Ypk6Fnj1DJxHJnnbtYNddYfz40Enyl4p8jnjgAejVyy/9FikUZnDbbXDffaGT5C+1a3LA11/7veIXLIC99059vEg+KSnxu6yOGuXblfJratfE3KBBcM45KvBSmKpX99sQ33tv6CT5SSP5wDZs8Bs1TZ0KhxwSOo1IGOvW+Z+D6dN174TyaCQfY08+6Vf+qcBLIdt+e7j6am11kAkayQfknJ8r3L+/v/uTSCH7/HM/2Fm2zM+4kV9oJB9Tr73mty7o2DF0EpHw9t4bzj4bhgwJnSS/aCQf0OmnwwUXQI8eoZOI5Ia5c+HMM/26kVq1QqfJHRrJx9CiRTB7Nlx8cegkIrnj8MP9hdcxY0InyR8q8oH07w/XXgu1a4dOIpJbbr7ZLw7Um/5oqF0TwNdfw4EHwocf+puDiMgvNm/2d0Z77DG/7YGoXRM7Q4dCly4q8CLlqVbN3/D7wQdDJ8kPGsln2aZNcMAB8MILcOSRodOI5KbvvvOLo+bNg/32C50mPI3kY2T8eGjUSAVeZFt22snfOGfQoNBJ4k8j+Szr0AFuuMFPnRSRrVu6FI47Dlas0J2jNJKPiblz/fzfLl1CJxHJfQcdBK1bw9NPh04SbyryWfTQQ3DddVCzZugkIvFw443+AmyBNwCqRO2aLNkybXLJEthjj9BpROJB0yk9tWtiYPhwOOssFXiRiqhWzb/7HTgwdJL4Smskb2adgP74XwrDnHP3lHNMEfAAUBNY7Zw7oZxjCnIkv3mz7y8+9RQcc0zoNCLxsmYNNG7stwKpXz90mjAyOpI3s2rAw8CpQAvgYjNrXuaYesBA4Ezn3KGA5o6UMnky7LILtGkTOolI/Oy8s5+NNnRo6CTxlE67pg2w1Dm3wjm3ERgNdC5zTDdgnHPuMwDn3FfRxoy3gQPh+uv9TYtFpOKuvx4eecQvJpSKSafINwBWlnr8afLPSmsK7Gpm08xslpl1jypg3H38McycCV27hk4iEl+HHeYXEU6YEDpJ/NSI8HmOBE4EdgBmmNkM59xHZQ8sLi7++fOioiKKiooiipCbBg2Cyy/XYg6Rqrr+ev+u+NxzQyfJvEQiQSKRiOS5Ul54NbNjgWLnXKfk476AK33x1cz6ALWdc39IPh4KvOycG1fmuQrqwuuGDdCwIcyYAU2ahE4jEm8//eT3sXnzTWjWLHSa7Mr0FMpZwIFm1sjMagFdgbJvmsYD7cysupltDxwDLKpMoHwybhwccYQKvEgUatWCK67Q7QErKmWRd86VADcAU4AFwGjn3CIzu9rMeiWPWQy8ArwPvAMMcc4tzFzseBg8GK65JnQKkfzRqxeMHAnr14dOEh9a8Zoh8+fDaafB8uVQI6orHyLC6af7iQyXXRY6SfZoxWsOeuQR6NlTBV4katdc498lS3o0ks+AtWv9Bdf334d99w2dRiS/bNrkV8BOmuSnVhYCjeRzzOjRcPzxKvAimVCjBlx1lX+3LKlpJJ8BbdpAcbHvHYpI9D79FFq1gpUrYYcdQqfJPI3kc8jcufDFF3DqqaGTiOSvfff1Ww8/+2zoJLlPRT5ijz7qL7hWrx46iUh+69VLc+bToXZNhNat8yvy5s1TP14k0zZtgv33h5dfhpYtQ6fJLLVrcsSzz0LbtirwItlQowZceaV/9yxbp5F8hI47Dvr0gbPPDp1EpDCsWAFHHeUvwObzJoAayeeABQv86lbNqBHJnkaN/Gy2sWNDJ8ldKvIRGTbMbymsFa4i2XXllf7nT8qndk0EfvrJ9+G1pbBI9m3Zgvjtt+HAA0OnyQy1awKbMAFatFCBFwmhVi249FJ47LHQSXKTinwEhg71bxlFJIwrr4QRI3QP2PKoyFfRJ5/ArFlw3nmhk4gUrkMO+WXOvPwnFfkqGjHC722dz9O3ROJAF2DLpwuvVbB5s+/DjxsHRx4ZOo1IYfv+e7/F96JFUL9+6DTR0oXXQBIJqFdPBV4kF+y4I5xzDjzxROgkuUVFvgqGD4cePUKnEJEtevTwP5d51DCoMrVrKum77/xbw6VLYY89QqcREfDF/aCD4OmnoXXr0Gmio3ZNAM8+CyeeqAIvkkvM/Mrz4cNDJ8kdGslXUrt2cPvt2oxMJNd88gkccQR89hnUrh06TTQ0ks+yJUvgo4/gtNNCJxGRsho29JMhxo8PnSQ3pFXkzayTmS02syVm1mcbx7U2s41mdm50EXPPiBFwySVQs2boJCJSHrVsfpGyXWNm1YAlQEfg38AsoKtzbnE5x70KrAcec849V85zxb5dU1LitzcthLvRiMTVunV+08D586FBg9Bpqi7T7Zo2wFLn3Arn3EZgNNC5nON6A2OBLysTJC6mTYO99lKBF8ll22/vtxp58snQScJLp8g3AFaWevxp8s9+Zmb7AF2cc4OASv22iYuRI+Gyy0KnEJFULrsMHn9cc+ajusVFf6B0r36rhb64uPjnz4uKiigqKoooQuatXeu3Fb7vvtBJRCSV446D9eth9mx/i8A4SSQSJBKJSJ4rnZ78sUCxc65T8nFfwDnn7il1zMdbPgV2B34AejnnJpR5rlj35EeOhDFjYOLE0ElEJB39+sGaNTBgQOgkVVOVnnw6Rb468CH+wuvnwEzgYufcoq0cPxyYmI8XXk86Ca65Bs4/P3QSEUnHsmXwm9/4OfNxng2X0QuvzrkS4AZgCrAAGO2cW2RmV5tZr/L+SmWC5LqVK2HOHDjzzNBJRCRdTZpA06aFvc+8Vrym6e67YcUKGDw4dBIRqYghQ+CVV/yW4HGV0XZNlOJa5J3zd54ZOtRfzBGR+Fizxq9tWb4cdtkldJrK0bYGGTZnDvz4I7RtGzqJiFTUzjvDySfD2LGhk4ShIp+GJ57wd4O3vF4BIJK/Lr20cG8monZNCps2wX77+btANWsWOo2IVMaPP8I++/g5840ahU5TcWrXZNDUqb7Iq8CLxNd228EFF8BTT4VOkn0q8imMGgXdu4dOISJV1b27/3mOWTOhytSu2Ya1a/1OdkuWwJ57hk4jIlXhHBxwADz3nL+pSJyoXZMh48f7KZMq8CLxZ+YvwI4aFTpJdqnIb8NTT0G3bqFTiEhULrkEnnnG3xeiUKjIb8Xq1fD229C5vJ3zRSSWmjf394N4883QSbJHRX4rxoyB00+HunVDJxGRKHXrVlg3E1GR3wq1akTyU9eu/uLrjz+GTpIdKvLlWL4cFi+GU04JnUREorbvvtCqVeHsTKkiX47Ro/2e8bVqhU4iIpnQrVvhLIzSPPlytGoFDz8M7duHTiIimfD1137O/MqVsNNOodOkpnnyEfrgA/jmG2jXLnQSEcmU3XaDDh3ghRdCJ8k8FfkyRo+Giy6Cavo/I5LXLr7Y/7znO7VrSnEODjrIL5aI293dRaRi1q6FBg38fWB33z10mm1TuyYi//ynX/p85JGhk4hIptWtC506+emU+UxFvpTRo/0cWt0cRKQwdO2a/y0btWuSNm/2NxOYPBlatAidRkSyYcMG2HtvWLjQ/zdXqV0TgX/8w98LUgVepHDUru33pxozJnSSzFGRT9rSqhGRwpLvLZu0iryZdTKzxWa2xMz6lPP1bmY2L/kx3cxaRh81czZt8r/JL7oodBIRybaOHWHpUr+dST5KWeTNrBrwMHAq0AK42MyalznsY6C9c+4w4M/Ao1EHzaQ33vD3cT3wwNBJRCTbataEc8+FZ58NnSQz0hnJtwGWOudWOOc2AqOB/9hl3Tn3jnPu2+TDd4AG0cbMrGee0ShepJBddFFhF/kGwMpSjz9l20W8JxCb/d02boTnn/d3cheRwtS+vd/HZtmy0EmiVyPKJzOzE4AewFZ3fikuLv7586KiIoqKiqKMUGHTpvmNivbfP2gMEQmoRg047zx/ba5v39BpIJFIkEgkInmulPPkzexYoNg51yn5uC/gnHP3lDmuFTAO6OScK/f3YS7Ok+/ZEw4+GG69NXQSEQkpkYBbboHZs0Mn+bWqzJNPp8hXBz4EOgKfAzOBi51zi0od0xB4HejunHtnG8+VU0V+40a/AGL2bGjYMHQaEQmppMTvZfPWW34Pq1yS0cVQzrkS4AZgCrAAGO2cW2RmV5tZr+RhdwK7An83szlmNrMyYbLt9dehaVMVeBGB6tX9zYLybWFUQW9rcMUV/gYhN90UOomI5II334TevWHevNBJ/lNG2zVRyqUi/9NPvlUzb56/56OISEmJXzMzbRo0axY6zS+0d00lTJ0KzZurwIvIL6pX97Nsxo4NnSQ6BVvkx471/TcRkdLOPz+/inxBtms2boR99vE3CdFFVxEpbcssm+nTc2erE7VrKiiRgCZNVOBF5NeqV/d72YwbFzpJNAqyyKtVIyLbkk8tm4Jr12za5Fs1774LjRsHjSIiOWpLnZg5Mze2PFG7pgLefNO3aVTgRWRratSALl3yYzRfcEV+7FjtOCkiqV1wQX6sfi2odk1JiZ8X/9ZbuXPVXERy05a9rebM8QukQlK7Jk0zZsCee6rAi0hqNWvCWWf5+03EWUEV+eee81OjRETSce65vm7EWcG0a5zzF1snTYJDDw0SQURiZsMGqF8flizxXYBQ1K5Jw+zZUKsWtGgROomIxEXt2tCpE4wfHzpJ5RVMkR83zm88ZJX6XSgiheq88+K9+rUg2jXO+R0nn3gCWrfO+rcXkRhbu9bvZbNiBey8c5gMateksHAhrF8PRx8dOomIxE3dunDCCTBxYugklVMQRX7LrBq1akSkMuK8YVlBtGuOOAL694cOHbL+rUUkD3zzDTRqBJ9/DjvskP3vr3bNNvzrX/DZZ9CuXegkIhJXu+wCbdrAK6+ETlJxeV/kx4/3q9aqVw+dRETi7Jxz4IUXQqeouLwv8s8/718cEZGq6NIFXnzR72kTJ2kVeTPrZGaLzWyJmfXZyjEPmtlSM5trZodHG7NyVq+GuXPhpJNCJxGRuGvQwO979cYboZNUTMoib2bVgIeBU4EWwMVm1rzMMacBTZxzBwFXA4MzkLXCJkyAU07xq9ayIZFIZOcbBaLzi698PjfI3vmdc078NixLZyTfBljqnFvhnNsIjAY6lzmmMzASwDn3LlDPzPaKNGklZLtVox+keMvn88vnc4PsFvkXXoDNm7Py7SKRTpFvAKws9fjT5J9t65jPyjkmq77/3t8F6owzQqYQkXzSrBnUqwezZoVOkr68vfA6eTK0betfEBGRqMStZZNyMZSZHQsUO+c6JR/3BZxz7p5SxwwGpjnnnkk+Xgx0cM6tKvNcYe/iLSISU5VdDFUjjWNmAQeaWSPgc6ArcHGZYyYA1wPPJH8prCmqvqQGAAAERUlEQVRb4KsSUkREKidlkXfOlZjZDcAUfHtnmHNukZld7b/shjjnXjKz083sI+AHoEdmY4uISDqyuneNiIhkV0YvvJrZLmY2xcw+NLNXzKzcy6BmVs/MxpjZIjNbYGbHZDJXVNI9v+Sx1cxstplNyGbGqkjn/MxsXzObmnzd5pvZjSGypiuuC/vSler8zKybmc1Lfkw3s5YhclZWOq9f8rjWZrbRzGJ1V+c0/30WmdkcM/vAzKalfFLnXMY+gHuA25Of9wH+spXjRgA9kp/XAHbKZK5sn1/y6zcDTwATQueO8vyA+sDhyc/rAh8CzUNn38r5VAM+AhoBNYG5ZbMCpwEvJj8/BngndO6Iz+9YoF7y8075dn6ljnsdmAScGzp3xK9fPWAB0CD5ePdUz5vpKZSdgceTnz8OdCl7gJntBBzvnBsO4Jzb5Jz7LsO5opLy/MCPdoHTgaFZyhWVlOfnnPvCOTc3+flaYBGB10hsQ2wX9qUp5fk5595xzn2bfPgOuftalSed1w+gNzAW+DKb4SKQzvl1A8Y55z4DcM59lepJM13k93TJWTbOuS+A8u533hj4ysyGJ9sZQ8ysToZzRSWd8wN4ALgNiNsFkHTPDwAz2x84HHg348kqJ5YL+yognfMrrSfwckYTRSvl+ZnZPkAX59wgIG6z+dJ5/ZoCu5rZNDObZWbdUz1pOlMot8nMXgVKj3QMX8z+u5zDyytyNYAjgeudc++ZWX+gL9CvqtmiUNXzM7MzgFXOublmVkSO/cOL4PXb8jx18aOn3yVH9JLDzOwE/Cy4fLvTQn98a3GLnPp5i8CWenkisAMww8xmOOc+2tZfqBLn3Mlb+5qZrTKzvZxzq8ysPuW/ffoUWOmcey/5eCz/+SIFFcH5HQecbWanA3WAHc1spHPusgxFrpAIzg8zq4F/3UY558ZnKGoUPgMalnq8b/LPyh6zX4pjclU654eZtQKGAJ2cc99kKVsU0jm/o4HRZmbA7sBpZrbROReHCQ/pnN+nwFfOuQ3ABjN7EzgM38svV6bbNROAy5Of/xb4VQFItgNWmlnT5B91BBZmOFdU0jm/3zvnGjrnDsAvJJuaKwU+DSnPL+kxYKFzbkA2QlXBzwv7zKwW/vUo+8M/AbgMfl7tXe7CvhyV8vzMrCEwDujunFsWIGNVpDw/59wByY/G+IHHdTEp8JDev8/xQDszq25m2+MnByza5rNm+GrxrsBr+BkXU4Cdk3++NzCp1HGHJU9wLvAcyav/uf6R7vmVOr4D8Zpdk/L88O9USpKv3RxgNn6EGDz/Vs6pU/J8lgJ9k392NdCr1DEP40dG84AjQ2eO8vyAR4Gvk6/THGBm6MxRv36ljn2MGM2uSff8gP+Hn2HzPtA71XNqMZSISB7L210oRURERV5EJK+pyIuI5DEVeRGRPKYiLyKSx1TkRUTymIq8iEgeU5EXEclj/x+loXDjy2pAjQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fef39c1ee90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x,sol,'-')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "<br>\n",
    "<br>\n",
    "Now we can proceed to compare the analytical solution with the numerical solution. <br>\n",
    "Remember, to execute the next cell you need to run first the sample utility."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAt4AAAHzCAYAAAAXThoqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XdUFOcaBvDno1gjxRYVFRC7Ith7xK4JxN4AEaOJvUUTjV4DSoxeS6Lm2kuwG2NXNGJXYokVe0dQ0FhBQQkI3/3DEruU3f12dp/fOZzD7M4Mz/oy8jK8MyuklCAiIiIiIv2yUB2AiIiIiMgcsPEmIiIiIjIANt5ERERERAbAxpuIiIiIyADYeBMRERERGQAbbyIiIiIiA1DaeAsh5gkh/hZCnHjPOh5CiGNCiFNCiJ2GzEdEREREpCtC5X28hRB1AMQDWCilrPCW520B7APQREoZLYTIK6W8Y+icRERERESZpfSMt5QyDMD996ziDWCVlDL62fpsuomIiIhIk4x9xrskgNxCiJ1CiENCiM6qAxERERERZYSV6gAfYAWgEoAGAHIC2C+E2C+lvPT6ikIIdTMzRERERGQ2pJQiI9sZe+N9HcAdKWUigEQhxB4AbgDeaLwBQOW8OmVcYGAgAgMDVcegDGL9tI310y7WTttYP+0SIkM9NwDjGDURzz7eZh2AOkIISyFEDgDVAZw1WDIyiKtXr6qOQJnA+mkb66ddrJ22sX7mSekZbyHEUgAeAPIIIaIABADIAkBKKWdLKc8JIbYAOAEgBcBsKeUZZYGJiIiIiDJIaeMtpfROwzoTAUw0QBxSxN/fX3UEygTWT9tYP+1i7bSN9TNPSu/jrUtCCGkqr4WIiIiIjJMQIsMXVxrDjDeZuV27dqmOQJnA+mkb66dd5lw7JycnCCH4wQ+9fjg5Oen8e9fY72pCRERE9IrIyEjeyYz0ToiM373knfs0lW9cwVETIiIisyCEYONNeveu77Nnj3PUhIiIiIjIWLHxJuXMeU7RFLB+2sb6aRdrR6Q9bLyJiIiIiAyAM95ERESkKcY+412+fHlMnz4dn3zyieoolAn6mPFm401ERESaYuyNd2Z5eHjg4MGDsLa2hpQSQghs3boV1atXVx3NrPDiSjJJnFPUNtZP21g/7WLt3i4iIhK+vqNQv34AfH1HISIi0qDb64IQAtOnT8eDBw/w8OFDPHjwgE23iWDjTURERCYhIiISjRv/giVLhmDXrlFYsmQIGjf+Jc3Nc2a3f87Z2Rk7duxAamoqfvzxRxQvXhy2traoWrUqoqOjAQD79u1DtWrVYG9vj+rVq2P//v2v7ONdZ/Tft11wcDDKli0LGxsbFC9eHLNnz37x3O7du1GkSBFMmDABH3/8MRwcHLBu3Tps3rwZpUqVQt68eTF27Nh0vU7KACmlSXw8fSlERERk6t71M9/HJ1AC8RKQL33ESx+fwDTtN7PbP+fk5CS3b98ux48fLytUqCAvXrwopZTyxIkT8t69e/LevXvS3t5eLlmyRKakpMhly5ZJe3t7ee/ePSmllB4eHnLevHlv7PdD223atElGRERIKaXcs2ePzJEjhzx27JiUUspdu3ZJKysr+cMPP8gnT57IOXPmyHz58kkfHx+ZkJAgT58+LbNnzy6vXr2artdqyt71ffbs8Qz1qzzjTURERCYhOjoVQM7XHs2JmJhUg2z/unnz5mHMmDEoXrw4AMDV1RX29vYICQlByZIl4e3tDQsLC3Ts2BGlS5fGhg0bXmzbv39/5M6dG/b29qhSpQoAfHC75s2bv3ib87p166JJkybYu3fvi31myZIFw4cPh6WlJTp27Ig7d+5g4MCByJEjB8qWLYuyZcsiPDw8Q6+V0oaNNynHOUVtY/20jfXTLtbuTQ4OFgASXns0AYUKpa3dyez2zz2/KO/atWsoVqzYG8/HxMTA0dHxlcccHR1fjKEAwNSpU3Hv3j3cv38fhw8fTtN2mzdvRs2aNZEnTx7Y29tj8+bNuHPnzot18+TJ8+Jt0LNnzw4AyJ8//4vns2fPjvj4+HS9VkofNt5ERERkEoKC/OHiEoB/m+cEuLgEICjI3yDbv0wIgaJFi+Ly5ctvPFeoUCFcvXr1lceioqLg4ODw3n2+b7ukpCS0bdsW3377LW7fvo379++jefPmJn33Fy1i403KeXh4qI5AmcD6aRvrp12s3ZucnR2xdWs/+PhMRP36AfDxmYitW/vB2dnxwxvrYPvnnje73bp1w8iRI3Hp0iUAwMmTJ3H//n18+umnuHjxIpYvX46UlBT89ttvOHv2LLy8vN673/dtl5SUhKSkJOTNmxcWFhbYvHkzQkND05Wb9M9KdQAiIiIiXXF2dsTixQHKtgfwYpxj8ODBSEpKQpMmTXD37l2ULl0aa9asQaFChbBx40b0798fvXr1QvHixRESEgJ7e/tXtn9d7ty537vd1KlT0a5dOyQlJcHLywstWrRIU853LZPu8Q10SLldu3bxzI2GsX7axvpplznXztTfQIeMA99Ah4iIiIhIo3jGm4iIiDSFZ7zJEHjGm4iIiIhIo9h4k3K8F622sX7axvppF2tHpD1svImIiIiIDIAz3kRERKQpnPEmQ+CMNxERERGRRrHxJuU4p6htrJ+2sX7axdoRaQ8bbyIiIiKNGDVqFDp37pyhbcPCwlCmTJlMZ3B2dsaOHTsyvZ+XWVhY4MqVKxnaVlevyxDYeJNy5vrOa6aC9dM21k+7WDvj5+Hhgdy5cyM5OVmn+03rW7u/3szWqVMHZ8+e1WkWXUnP29Vr6XW9jo03ERERkY5FRkYiLCwMFhYWWL9+vZIM6WlmVUvPxbJael2vY+NNynFOUdtYP92IiIiEr+8o1K8fAF/fUYiIiDTI1zVE/VS9NlPHY8+4LVy4EDVr1oS/vz+Cg4Nfea5r167o27cvPD09YWNjg5o1ayIiIuLF8wMHDkTRokVha2uLqlWrIiws7K1fw9PTE9OmTXvlMTc3N6xbtw716tWDlBIVKlSAjY0Nfv/9d+zevRtFihR5se7169fRpk0b5M+fH/ny5UP//v0BAFeuXEHDhg2RN29e5M+fH76+vnjw4EGaXvemTZtQrlw52NjYoEiRIvjpp59ePDdnzhyUKFECefPmRcuWLXHjxo237qN+/fqYP3/+i+UFCxagbt26AJCm13Xu3DnUr18f9vb2cHV1xYYNG14896F/e31j401EpFhERCQaN/4FS5YMwa5do7BkyRA0bvyLSTSopvzaiN5n4cKF8PX1hbe3N7Zs2YLbt2+/8vxvv/2GUaNGITY2Fi4uLhgxYsSL56pVq4YTJ07g/v378Pb2Rrt27ZCUlPTG1+jSpQsWLVr0Yjk8PBwxMTHw9PTE7t27AQAnT57EgwcP0K5dOwD/ni1OTU2Fp6cnnJ2dERUVhejoaHTs2BHA07PPw4cPx82bN3H27Flcv34dgYGBaXrd3bt3x5w5c/DgwQOcOnUKDRo0AADs2LEDw4cPx8qVK3Hjxg0ULVr0xddLi+e5P/S6njx5Ai8vLzRr1gy3b9/G1KlT4ePjg4sXL77Y1/v+7fXNymBfiegdOKeobaxf5o0cGYzLl0cByPnskZy4fHkURo6ciMWLA9K8Hykl7jy6g8i4SETFReFm/E3EJcYh7p84xCbGIu6fOMQlxiFFpryy3djFY2EhLGCT1QZ2We1gm80WtlltYZfNDvlz5kdR26JwtHNE/pz5YSHSd75GV6+N3sRj7/3EqMyPI8iAjN0rPCwsDFFRUWjfvj3s7e1RvHhxLF26FAMGDHixTqtWrVC5cmUAgI+PDwYPHvziOW9v7xefDxo0CEFBQTh//jxcXV1f+Tqff/45evbsicuXL8PFxQWLFy9Ghw4dYGlp+e9reMcIx8GDB3Hjxg2MHz8eFhZPj+tatWoBAFxcXODi4gIAyJMnDwYNGoTRo0en6bVnyZIFp0+fhqurK2xtbeHu7g4AWLp0Kbp16wY3NzcAwNixY2Fvb4+oqCgULVo0Tft+2bte1/79+5GQkIChQ4cCeHr23NPTE8uWLcP3338P4P3/9vrGxpuISLHo6FT825g+lxMxMalvXf9x8mOcuX0GJ/4+gZO3TuL07dO4GnsV1+KuIbt1djjaOqKobVEU+KgA7LLZwTarLRxtHWGbzRY2WW1gbWH9xj5TZAoe/PPgaYP+rFmPfhiNHVd3ICouCpGxkXjwzwMUtikMRztHlM1bFhU+rgDXj11RPn95fJTlI528NiJdyWjTrAsLFy5EkyZNYG9vDwDo1KkTFixY8ErjXaBAgRef58iRA/Hx8S+WJ06ciPnz578YxXj48CHu3LnzxtfJmjUrOnTogMWLF+P777/HsmXLsGrVqjRlvH79OhwdHV803S+7desWBgwYgL179yI+Ph4pKSnInTt3mva7atUqBAUFYejQoXBzc8O4ceNQvXp1xMTEvGh2ASBnzpzIkycPoqOjM9R4v8uNGzdeGTsBAEdHR0RHR79Yft+/vb6x8Sbldu3axTM3Gsb6ZZ6DgwWABLzaoCagUCELpMpUnLp1Cnsj9yLsWhiO3TiGyLhIlMxTEq75XVHh4wpoVKwRnO2cUcS2yDsb4HdJT/0eJz/GtQfXcDX2Kk7fOo191/dh1pFZOHP7DArmKoiKBSqidpHaqOtYF+4F3GFlYfXe10aZw2PPOCUmJmLFihVITU1FwYIFAQBJSUmIjY3FyZMn3zhr/bq9e/diwoQJ2LlzJ8qWLQsAyJ079zvP8Pr5+aFz586oXbs2cubMierVq6cpZ5EiRRAVFYXU1NQ3mu/hw4fDwsICp0+fhq2tLdatW4d+/fqlab+VK1fG2rVrkZKSgl9++QXt2rVDVFQUChUqhMjIf0fMEhIScPfuXRQuXPiNfeTMmROPHj16sXzz5s00fW0AKFSoEK5du/bKY1FRUShVqlSa96FPbLyJiBQLCvLHgQMBz0YycgAF9iN3tWGI8bBG3vF5kTdHXtQtWhdNijXBiLojUDJPSWSxzGLwnNmts6NknpIomackmrg0efH4k9QnuHTvEo7EHEFYVBjmHZuHqLgo1ChcA+XbucLh0peI/msmIG0AJMDFJQBBQWn7IU6kNWvWrIGVlRXCw8Nhbf3vX5fatWuHhQsXYsKECe/dPj4+HtbW1siTJw+SkpIwbtw4PHz48J3r16hRAxYWFhg8ePAb9/cuUKAArly5gmLFir2xXbVq1VCwYEEMGzYMgYGBsLS0xJEjR1CrVi08fPgQdnZ2yJUrF6Kjoz+Y+bnk5GT8/vvvLy5czJUr14uxl06dOsHb2xve3t4oVaoUhg8fjho1arxxdhoA3N3dsXr1anTr1g3R0dGYN2/eK2ep3/e6qlevjhw5cmD8+PH4+uuvERYWho0bN6Z5Rl3feMqBlOMZG21j/TIvv0NeDJ1fGsUH1UOW72yR64vP4eXtgt61euFMnzO40O8C5rWYh64Vu6J8/vI6bbp1UT8rCyuUzlsaPhV8MMNzBk71PoWIARHoW60vLLNb4KNOx5B1REEU6OmGOj26YlVIFzg7O2Y+vJnjsWecFi5ciC+++AIODg7Inz//i4++fftiyZIlSE19/5hV06ZN0bRpU5QsWRLOzs7IkSPHW5vTl/n5+eHUqVPw9fV95fHAwED4+fkhd+7cWLly5SvPWVhYYMOGDbh48SKKFi2KIkWKYMWKFQCAgIAAHDlyBHZ2dvDy8kKbNm1e2fZ9t/NbtGgRnJ2dYWdnh9mzZ2Pp0qUAgIYNGyIoKAitW7eGg4MDIiIisHz58rfuc9CgQbC2tkaBAgXQtWvXdL0ua2trbNiwAZs2bULevHnRt29fLFq0CCVKlPhgdkMQ6blvojETQkhTeS1EZPrik+Kx9txaLDu1DHsj96KqQ1V4lvDEZyU/Q8k8JVXH07nI2EiEXAzBxgsbERYVhqoOVdGxXEe0LdsW9tntVccjjRFCpOu+z6Zu0aJFmDNnDvbs2aM6ikl51/fZs8cz1MGz8SblOKeobaxf2iWnJCP0ciiWnFyCTRc3oU7ROvB29cZnJT6DbTZbJZlU1C8hKQGhl0Ox9NRShF4ORQPnBvBx9YFnSU9ks8pm0CxaZs7HHhvvfz169AgNGzZE37594ePjozqOSdFH480ZbyIiPTt/5zxmHJ6BpSeXokSeEvBx9cHU5lORN0de1dGUyJklJ1qVaYVWZVohLjEOq8+uxszDM/HVhq/Qpkwb9K7aGxULVlQdk8johYaGonXr1mjSpAk6deqkOg6lAc94ExHpwZPUJ1h/fj2mH5qOk7dOonvF7uhWqRuK2b95MRA9Ff0gGsHHgzHryCw42DigT9U+aFe2HbJaZVUdjYwMz3iTIXDU5D3YeBORMbjz6A5mHp6JWUdmwdHWEb2r9kabMm3YPKbDk9QnCLkQgmmHpiH873B84f4F+lTrg8I2b952jMwTG28yBH003ryrCSm3a9cu1REoE1i/p64/uI5BfwxCyV9K4mrsVWzotAFhX4TB29XbqJtuY6yflYUVWpRugdDOodjbdS8eP3mMCjMq4Mv1X+Li3Ysf3oGZMMbaEdH7sfEmIsqEi3cvovv67qgwowIsLSxxqvcpzP18LtwLuKuOZhJK5imJyc0m42K/i3CwcUCt+bXQcWVHhN8MVx2NiCjdOGpCRJQBF+5eQMCuAGy7sg19q/ZF32p9kSdHHtWxTN7Dfx5i9pHZ+OnAT6hYoCKC6gfxQkwz5OTk9Mq7IBLpg6OjI65evfrG45zxBhtvIjKMmIcxGLVrFFafW42va3yNftX7pftt2inzEp8kYt7Refhh7w+o71QfQfWD4JLbRXUsIjIDnPEmTeOcoraZS/1iE2Px3bbv4DrDFbbZbHG+73l8V/c7zTfdWq1fNqts6FOtDy72u4iy+cqi+tzq6BPSBzfjb6qOZjBarR09xfqZJzbeRETvkZySjJ/3/4ySv5TEnUd3EN4zHOMbj0fu7LlVRyMAH2X5CP/55D841/ccslplRbnp5RC4KxCPkh+pjkZE9AaloyZCiHkAPAH8LaWs8J71qgLYB6CDlHL1O9bhqAkR6dSOiB3ot7kfCtsUxuSmk1EmXxnVkegDImMj8e22b3Hw+kH83PRntCzdEkJk6C/CRERvpdkZbyFEHQDxABa+q/EWQlgA2ArgMYD5bLyJSN+uxV3D4NDBOBRzCD83/RktSrVg86YxL//SNLXZVJTKW0p1JCIyEZqd8ZZShgG4/4HV+gFYCeCW/hORCpxz0zZTql9SShJ+3Psj3Ge5o2y+sjjd+7TJnzE1pfq9rIFzAxzvcRxNXZqi9vzaGLp1KBKSElTH0ilTrZ25YP3Mk1HPeAshCgFoKaWcAcB0f/IRkXKHog+h8uzK2HdtHw59eQiBHoHIYZ1DdSzKBGtLa3xd82uc7HUS0Q+j4TrDFduvbFcdi4jMmPLbCQohHAFseNuoiRBiBYCJUsq/hBC/AtgopVz1jv3ILl26wMnJCQBgZ2cHd3d3eHh4APj3N0suc5nLXH55+XHyY/hP9kfolVBM6z0Nncp3wu7du40mH5d1t/zI4RF6buwJ10eu6FW1FzybeBpVPi5zmcvGuXz8+HHExsYCAK5evYoFCxZoc8Yb+GDjfeX5pwDyAkgA8JWUcv1b1uWMNxGly97Ivei2vhsqFayEqc2nIn/O/KojkZ49+OcBvt36LUIuhmDGZzPgWdJTdSQi0hjNzng/I/COMRIpZbFnH854Oufd+21NN2nb898uSZu0WL9HyY/Qb1M/dFzVEeMbj8fytsvNtunWYv0ywyarDWZ6zsTClgsx8I+B8F3ti9jEWNWxMsTcamdqWD/zpLTxFkIsxdPbBJYUQkQJIboKIXoIIb56y+o8nU1EmXb85nFUmV0Fdx/fxalep9CydEvVkUiB+s71caLXCdhls4PbTDfsidyjOhIRmQHloya6wlETInqfVJmKyQcmY2zYWPzc9Gf4VvBVHYmMRMiFEHTf0B3dKnZDQL0AWFtaq45EREZMs/fx1iU23kT0LjEPY+C/1h/xSfFY0noJnO2dVUciI/N3/N/ouq4r7j6+iyWtl6B47uKqIxGRkdL6jDeZOc65aZux1y/kQggqzaqE2kVqY0/XPWy6X2Ps9TOUjz/6GCHeIfBx9UHNeTWxKHyR6kgfxNppG+tnnqxUByAi0oeU1BR8v/N7LDyxECvbr0SdonVURyIjJ4RA/+r94eHkgbYr2uLPa39icrPJyGaVTXU0IjIRHDUhIpNzK+EWvFd5Q0JiWZtlZnvHEsq4B/88QNd1XREVF4Xf2/0OJzsn1ZGIyEhw1ISI6Jn91/aj8uzKqO5QHaG+oWy6KUNsstpgZbuV6FS+E6rPrY7NFzerjkREJoCNNynHOTdtM1T9IiIi4es7CvXrB8DXdxQiIiJfeV5KiV8O/oIWy1tg+qfTMabhGFhaWBokm5bx+Hs3IQS+rvk1VrZbiS83fImAnQFIlalvrPeh7019Ye20jfUzT5zxJiKjFxERicaNf8Hly6MA5ASQgAMHArB1az84Ozvinyf/oFdILxy5cQQHuh9AMftiqiOTCanrWBeHvzqM9r+3x/HfjmNxq8XIlTUXgA9/bxIRvYwz3kRk9Hx9R2HJkiF42tg8lwAfn4n4aVYvtFnRBvly5MPCVgvxUZaPVMUkE5eUkoQ+IX1wMPog1ndaDyc7p/d+by5eHKAqKhHpEWe8icikRUen4tXGBgBy4kLcDVSfWx31HOthZfuVbLpJr7JYZsFsr9noVrEbas6ribCosHd+b8bEvDmSQkTExpuU45ybthmifg4OFgASXn2w1AqccF+EMQ3G4IcGP8BC8L+zjODxlz5CCAyoMQDBLYLR+rfWSCwdjje+N5GAQoX0//3I2mkb62ee+JOKiIxeUJA/XFwC8LTBkUCtMbD8vBuWf74U3q7eitOROWpavCn2dN2DGyXCYdvOAxAPnz2TABeXAAQF+asLR0RGizPeRKQJERGRGDFyPvbk2IyEvJEI8V6HWuVrqI5FZu7e43toFtwcty4nwOlYCxQumAVBQf68sJLIhGVmxpuNNxFpwuPkx/Bd44t7j+9hbYe1sM1mqzoSEQAg8UkifFb74N7je1jTYQ3sstmpjkREesSLK0nTOOembYao373H99BkcRNYW1jjD58/2HTrEI+/zMtmlQ0r2q5A+Xzl8cmvnyD6QbRBvi5rp22sn3li401ERi0qLgp15tdBtULVsLTNUmS1yqo6EtEbLC0sMbX5VPhW8EWt+bVw+tZp1ZGIyAhx1ISIjNapW6fQfElzDKoxCF/X/Fp1HKI0WRS+CEO2DsGq9qtQp2gd1XGISMc44w023kSm5lD0IXgt88JPTX/inUtIc7Zc2gLfNb5Y2nopGrs0Vh2HiHSIM96kaZxz0zZ91G9P5B58uvRTzPKcxaZbz3j86UfT4k2xuv1q+Kz2wdpza/XyNVg7bWP9zJOV6gBERC/jmUIyFXUd62KTzyZ4LvXEo+RH/CWSiDhqQkTGY83ZNeixsQfWdFiD2kVrq45DpBOnbp1C08VNEVgvEF9W/lJ1HCLKpMyMmvCMNxEZhSUnlmBw6GBs9tmMyoUqq45DpDPl85fHri670HhRYzxMesgLhYnMGGe8STnOuWmbLuq3MHwhvt32Lbb7bWfTbWA8/gyjRJ4S2NN1D2YcnoEJf07QyT5ZO21j/cwTz3gTkVKLTyzGsG3DsN1vO8rkK6M6DpHeFLUtip1ddsIj2AOWFpY8801khjjjTUTKLD25FENCh2Cb3zaUzVdWdRwig7gWdw0eCzzQr1o/DKwxUHUcIkonzngTkeb8duo3DA4djK2dt7LpJrNSxLbIv2e+hSX6Ve+nOhIRGQhnvEk5zrlpW0bq9/vp3zFwy0CE+oaifP7yug9FacbjT43nYyc/HfgJ0/6alqF9sHbaxvqZJ57xJiKDWn12Nfr/0R9bfLfA9WNX1XGIlHG0c3xl5rtnlZ6qIxGRnnHGm4gMZsulLfBb64ctvlvgXsBddRwio3Dl/hXUC66HcQ3HwaeCj+o4RPQBnPEmIqO379o+dF7TGWs6rGHTTfSSYvbFsMV3CxosaACbrDbwKuWlOhIR6QlnvEk5zrlpW1rqF34zHK1+a4VFrRbxHSmNDI8/41A2X1ls6LQB3dZ3w86InWnahrXTNtbPPLHxJiK9unD3ApovaY5pn05D0+JNVcchMlpVHapiRbsV6LCyAw5FH1Idh4j0gDPeRKQ31+Kuoe6vdfF9ve/xRcUvVMch0oQN5zfgyw1fYrvfdpTLX051HCJ6TWZmvHnGm4j04nbCbTRe1Bj9q/dn002UDl6lvDCpySQ0XdwUEfcjVMchIh1i403Kcc5N295Wv0fJj+C1zAttyrTh22IbOR5/xsmngg+G1RmGZkua4c6jO29dh7XTNtbPPLHxJiKdepL6BB1XdkTpvKXxQ4MfVMch0qy+1fqiVelW+HzZ53iU/Eh1HCLSAc54E5HOSCnRc2NPXI27io2dNsLa0lp1JCJNS5Wp8Fvjh4TkBKxstxKWFpaqIxGZPc54E5FRGLN3DA7FHMLKdivZdBPpgIWwwPwW8/Hwn4fov7k/eIKJSNvYeJNynHPTtuf1Cz4ejHnH5iHEOwS5suZSG4rSjMef8ctimQWr2q9C2LUw/PfP/754nLXTNtbPPLHxJqJM23JpC4ZtG4bNPptRMFdB1XGITI5tNlts9tmMmYdnYlH4ItVxiCiDOONNRJly8u+TaLiwIdZ0WMN3pSTSszO3z8Aj2AO/t/sd9ZzqqY5DZJY4401ESvwd/ze8lnlhSrMpbLqJDKBsvrJY2mYpOqzsgEv3LqmOQ0TpxMablOOcmzYlPklEy99aop6sh06unVTHoQzi8ac9jYo1QqBHIBqMaoD7j++rjkMZxGPPPLHxJqJ0k1Lii3VfwNHWEf7u/qrjEJmdnlV6oppDNbRf2R7JKcmq4xBRGnHGm4jSLWh3EDZe3IhdXXYhu3V21XGIzFJKagq8lnnB0dYR0z+bDiEyNHJKROnEGW8iMpgVp1dg7rG5WNdxHZtuIoUsLSyxvO1y7I3ai1/++kV1HCJKAzbepBzn3LRR3l5+AAAgAElEQVTjUPQh9NnUB+s7rkeBjwoAYP20jvXTrl27dsEmqw02em/E2LCx2Hxxs+pIlA489swTG28iSpOb8TfRekVrzPWaC7cCbqrjENEzTnZOWNluJbqs7YILdy+ojkNE76F0xlsIMQ+AJ4C/pZQV3vK8N4ChzxYfAuglpTz5jn1xxptIT5JSktBgQQM0LtYYAR4BquMQ0VvMPjIbkw9MxoHuB2CT1UZ1HCKTlZkZb9WNdx0A8QAWvqPxrgHgrJQyTgjRDECglLLGO/bFxptIT3qH9EbMwxis7rAaFoJ/KCMyVj029MCtR7ewqv0qHqtEeqLZiyullGEA3nkTUinlASll3LPFAwAcDBKMDIpzbsZt7tG52Hl1Jxa2WvjWH+Ssn7axftr1ttpNbT4Vf8f/jTF7xhg+EKULjz3zpKVfh7sD4JUjRAZ04PoBDN8+HGs7rOWfrok0IKtVVqxqvwqzjszCxgsbVcchotcov4+3EMIRwIa3jZq8tE59AP8DUEdK+dYz5EII2aVLFzg5OQEA7Ozs4O7uDg8PDwD//mbJZS5zOW3Ldx/dxYBzAzDTcyY+ivlIeR4uc5nLaV8+c/sMAq8GYm/Xvbhx6obyPFzmspaXjx8/jtjYWADA1atXsWDBAm3OeAMfbryFEBUArALQTEp5+T374Yw3kY4kpSTBI9gDzYs3x8h6I1XHIaIMmHt0Libtn4SD3Q/yL1ZEOqTZGe9nxLOPN58QoiieNt2d39d0k7Y9/+2SjMeQ0CHInzM/Rnwy4oPrsn7axvpp14dq171Sd9RzrIdu67uBJ6aMD48986S08RZCLAWwD0BJIUSUEKKrEKKHEOKrZ6uMBJAbwHQhxDEhxF/KwhKZid9O/YZNFzchuGUw74pApHGTm01GxP0ITD04VXUUIoIRjJroCkdNiDLv3J1zqPtrXYT6hqJiwYqq4xCRDkTcj0CNeTWwtsNa1CxSU3UcIs3T+qgJERmBhKQEtF3RFmMbjmXTTWRCnO2dMddrLjqs7IDbCbdVxyEya2y8STnOuaknpUTPkJ6oXKgyulXslq5tWT9tY/20Kz218yrlBW9Xb/iu8UVKaor+QlGa8dgzT2y8iQhzjs7B8ZvHMf3T6RAiQ389IyIj90ODH5D4JBE/7PlBdRQis8UZbyIzdyTmCJotaYawrmEolbeU6jhEpEc3Ht5A5dmVEdwyGE1cmqiOQ6RJnPEmogyJS4xD+5XtMe3TaWy6icxAwVwFsbTNUvit8UPMwxjVcYjMDhtvUo5zbmo8n+tuUqwJ2pdrn+H9sH7axvppV0Zr5+Hkgd5Ve8N3Nee9VeKxZ57YeBOZqV+P/4pTt07hp6Y/qY5CRAY2ou4IpMpUjAsbpzoKkVnhjDeRGTp7+yw+Cf4Eu7rsQrn85VTHISIFrj+4jiqzq2BV+1WoXbS26jhEmsEZbyJKs8Qniei4qiN+bPAjm24iM1bYpjDmeM2B92pv3H98X3UcIrPAxpuU45ybYQ0JHYJSeUqhe6XuOtkf66dtrJ926aJ2XqW80LJUS3Tf0B38q7Fh8dgzT2y8iczImrNrEHIxBLO9ZvN+3UQEABjfeDwi7kdg5uGZqqMQmTzOeBOZiWtx11BlThWs67gONQrXUB2HiIzIhbsXUHt+bWz3244KH1dQHYfIqHHGm4jeKyU1BZ3XdMbA6gPZdBPRG0rmKYkJjSfAZ7UPEp8kqo5DZLLYeJNynHPTv5/2/4RUmYpva3+r832zftrG+mmXrmvXxa0LSuctjeHbh+t0v/R2PPbMExtvIhN3/OZxTNg3AQtbLYSlhaXqOERkpIQQmPnZTKw4vQLbrmxTHYfIJHHGm8hEREREYuTIYERHp8LBwQJBQf4oUDg/qs6piqG1h6KzW2fVEYlIA7Ze3oov1n+B8J7hyJ0991v/b3F2dlQdk0iZzMx4s/EmMgEREZFo3PgXXL48CkBOAAlwcQlA3R/v4ZFFApa3Wc67mBBRmg38YyBiHsZgXOXxaNLkf2/837J1az8232S2eHElaRrn3DJv5Mjgl34wAkBOXJb1sOLEKsz4bIZem27WT9tYP+3SZ+3GNhyLM7fPwHfC12/+33J5FEaODNbb1zYXPPbMExtvIhMQHZ2Kf38wAsh+F2jZCyXOfI7c2XMry0VE2pTdOjuWtF6Cw7m3AHa3X3s2J2JiUpXkItI6Nt6knIeHh+oImufgYAEg4dmSBDx7AqdboXyO4nr/2qyftrF+2qXv2rkVcEP5uNpAKx9ApLz0TAIKFWL7kFk89swTjxwiExAU5A8XlwAACYDrMiDfaThfsURQkL/aYESkaSsGzUS2rFFAjfHPHnk6483/W4gyho03Kcc5t8xzdnbE1q390Mrve1h7fYVmiXWxfcsgg1z8xPppG+unXYaoXXGXYvij11JkbfgDqn3WGz4+E3lhpY7w2DNPVqoDEJFuODkVRXKzCxhW4GuMrj9adRwiMhH1KtTFz/9MRLBjMIK/+BNWFmwdiDKKtxMkMhHBx4Mx+cBk/PXlX8himUV1HCIyIakyFU0WNUFD54b4ru53quMQKcX7eIONN5m36w+uo+KsitjWeRvcCripjkNEJigyNhJV5lTBDr8dcP3YVXUcImV4H2/SNM65ZY6UEt3Xd0f/av2VNN2sn7axftpl6No52jlibMOx8F/nj+SUZIN+bVPEY888sfEm0ri5R+fizqM7GFZnmOooRGTiulXshvw582Ns2FjVUYg0iaMmRBr2/E+/O7vsRPn85VXHISIz8Hy0bWvnrXAv4K46DpHBcdSEyAxJKdFtfTcMrjmYTTcRGUxhm8KY2HgiuqztwpETonRi403Kcc4tY349/itiE2MxpNYQpTlYP21j/bRLZe383PzgkMsB4/8c/+GV6a147JknNt5EGhTzMAZDtw3FvM/n8Z66RGRwQgjM9JyJnw/8jLO3z6qOQ6QZnPEm0qDWv7VG2Xxl8UODH1RHISIzNu2vaVh6ain2dt0LC8FzeWQeOONNZEZWnVmFs3fO4j+f/Ed1FCIyc72q9oKAwPRD01VHIdIENt6kHOfc0u7e43vot7kf5nrNRTarbKrjAGD9tI710y5jqJ2FsMC8z+chcFcgImMjVcfRFGOoHxkeG28iDRkcOhhty7ZF7aK1VUchIgIAlMpbCoNrDkaPjT3AkU+i9+OMN5FGhF4OxVcbvsKp3qfwUZaPVMchInohOSUZ1eZWw6Aag+Dn5qc6DpFeccabyMTFJ8Wjx8YemOU5i003ERkda0trzPt8Hr7Z+g3+jv9bdRwio8XGm5TjnNuHBewMQN2iddG0eFPVUd7A+mkb66ddxla7SgUrwd/NH4O2DFIdRROMrX5kGGy8iYzcsRvHsPjkYkxqMkl1FCKi9wrwCMCB6wcQejlUdRQio8QZbyIjlpKagprzaqJXlV7oWrGr6jhERB+0+eJm9N3cF6d6nUJ26+yq4xDpHGe8iUzUjMMzkN06O/zd/VVHISJKk+YlmqNKoSr4YQ/f4IvodWy8STnOub1d9INoBO4KxMzPZkKIDP1ibRCsn7axftplzLWb3HQyZh+djdO3TquOYrSMuX6kP2y8iYzUgD8GoHfV3iiTr4zqKERE6VIwV0GM9hiNHht7IFWmqo5DZDQ4401khDZe2IhBWwbhZK+TRvMOlURE6ZEqU1FrXi10r9Qd3St1Vx2HSGc4401kQhKSEtB3U1/M+GwGm24i0iwLYYFZnrMwfPtw3tub6Bk23qQc59xeFbgrEHUd66JRsUaqo6QJ66dtrJ92aaF2bgXc4O/uj8Ghg1VHMTpaqB/pHhtvIiNy6tYpLAhfwHt2E5HJCKgXgLCoMOyM2Kk6CpFySme8hRDzAHgC+FtKWeEd60wF0BxAAgB/KeXxd6zHGW/SNCklPBZ4oGO5juhVtZfqOEREOrP23FoM3z4c4T3DYW1prToOUaZoecb7VwDvfA9sIURzAC5SyhIAegCYaahgRIa29ORSxCfF46vKX6mOQkSkUy1KtYCTnROmHJyiOgqRUkobbyllGID771mlBYCFz9Y9CMBWCPGxIbKR4XDODYhLjMM3W7/B9E+nw9LCUnWcdGH9tI310y4t1U4IganNp2Jc2Dhcf3BddRyjoKX6ke6oPuP9IQ4Arr20HP3sMSKTErgrEJ+V+AzVC1dXHYWISC+K5y6O3lV7Y0joENVRiJSxUh1Al/z9/eHk5AQAsLOzg7u7Ozw8PAD8+5sll41v2cPDw6jyGHr5xN8nELwuGMEtgvGcMeX70LK510/ry6wflw25PKzOMBQbVAyTkiZhsPdg5Xm4zOW0LB8/fhyxsbEAgKtXryIzlL+BjhDCEcCGt11cKYSYCWCnlPK3Z8vnANSTUr5xQ1BeXElaJKXEJ8GfwNfVFz2q9FAdh4hI79afX4+h24YivGc4slhmUR2HKN20fHElAIhnH2+zHoAfAAghagCIfVvTTdr2/LdLc7T4xGIkPknU9Lu6mXP9TAHrp11arZ1XSS+42Ltg8oHJqqMopdX6UeYoHTURQiwF4AEgjxAiCkAAgCwApJRytpRykxDiUyHEJTy9nWBXdWmJdCs2MRZDtw3Fuo7rNHdBJRFRRgkhMKXZFFSfWx2dyndCEdsiqiMRGYzyURNd4agJac3APwbiUfIjzPaarToKEZHBBe4KxJnbZ7Ci3QrVUYjSJTOjJmy8iRQ4fes06i+ojzN9ziBvjryq4xARGdzj5McoM60MglsGw8PJQ3UcojTT+ow3mTlzm3OTUmLQlkH4zyf/MYmm29zqZ2pYP+3Seu2yW2fHhMYTMOCPAUhJTVEdx+C0Xj/KGDbeRAa2/vx6RD+MRq8qfFt4IjJvbcu2hX02e8w5Okd1FCKD4KgJkQElPklEuenlMMtzFhoVa6Q6DhGRcuE3w9FkcROc63MO9tntVcch+iCOmhBpxM/7f4Zrflc23UREz7gVcEPr0q0RsCtAdRQivWPjTcqZy5xbzMMYTNo/CZOaTFIdRafMpX6mivXTLlOqXVCDICw/tRynb51WHcVgTKl+lHZsvIkMZNi2Yfiq8ldwye2iOgoRkVHJmyMvRn4yEgO3DATHRsmUccabyAAOXD+ANiva4Fyfc8iVNZfqOERERic5JRnus9wxpsEYtCzdUnUconfijDeREUuVqei/uT/GNRzHppuI6B2sLa0xpdkUDA4djMQniarjEOkFG29SztTn3BafWAwLYQGfCj6qo+iFqdfP1LF+2mWKtWtUrBFc87ti8oHJqqPonSnWjz6MjTeRHj1KfoQRO0bg56Y/w0LwcCMi+pAJjSdg4r6JuJVwS3UUIp3jjDeRHgXtDsLp26exvO1y1VGIiDRj0B+D8E/KP5j+2XTVUYjekJkZbzbeRHoS8zAGrjNccfjLw3C2d1Ydh4hIM+49vofS/yuNnV12olz+cqrjEL2CF1eSppnqnNvIHSPRvWJ3k2+6TbV+5oL10y5Trl3u7LkxvO5wfLP1G9VR9MaU60fvxsabSA/Cb4Yj5GIIhtcdrjoKEZEm9a7aGxfvXUTo5VDVUYh0hqMmRDompUSjRY3Qtkxb9KraS3UcIiLNWntuLUbuHInjPY7D0sJSdRwiABw1ITIqIRdDcOPhDXxZ+UvVUYiINK1FqRbInT035h+brzoKkU6w8SblTGnOLTklGd9s/QYTm0yElYWV6jgGYUr1M0esn3aZQ+2EEPipyU8I2BWAh/88VB1Hp8yhfvQmNt5EOjT7yGwUsSmC5sWbq45CRGQSKheqjMYujfHfP/+rOgpRpnHGm0hH4hLjUPJ/JbG181ZU+LiC6jhERCbj+oPrcJvphuM9jqOIbRHVccjMccabyAiM/3M8Pi3xKZtuIiIdK2xTGD0q90DArgDVUYgyhY03KWcKc27RD6Ix88hMjPYYrTqKwZlC/cwZ66dd5la7obWHIuRiCE7dOqU6ik6YW/3oKTbeRDoQuCsQ3St2559AiYj0xDabLb6r8x2GbRumOgpRhnHGmyiTztw+A49gD5zvex722e1VxyEiMln/PPkHZaaVwa8tfkU9p3qq45CZ4ow3kULfbf8OQ2sPZdNNRKRnWa2y4ocGP+Dbbd+CJ9tIi9h4k3JannMLiwpD+M1w9KnWR3UUZbRcP2L9tMxca9exfEc8SX2ClWdWqo6SKeZaP3PHxpsog6SU+GbrNwiqH4RsVtlUxyEiMgsWwgL/bfRfDN8xHEkpSarjEKULZ7yJMmj12dUYvXs0jvY4CgvB32GJiAyp6eKm8Crphb7V+qqOQmYmMzPebLyJMiA5JRnlZ5TH1GZT0bR4U9VxiIjMzvGbx9FscTNc6HcBNlltVMchM8KLK0nTtDjnNv/YfBSxKYImLk1UR1FOi/Wjf7F+2mXutXMv4I7GLo0xcd9E1VEyxNzrZ67YeBOl06PkRxi9ZzTGNRoHITL0Cy8REelAUP0gTDs0DbcSbqmOQpQmHDUhSqcJf07AweiDWNle21fUExGZgn6b+sHKwgo/N/tZdRQyE5zxBhtvMoy4xDiU+KUEdvvvRpl8ZVTHISIyezfjb6Lc9HI41uMYitoWVR2HzABnvEnTtDTnNnHfRHxW8jM23S/RUv3oTayfdrF2TxX4qAB6VO6B0btHq46SLqyfeWLjTZRGf8f/jemHpyOwXqDqKERE9JJvan2DdefX4fyd86qjEL0XR02I0mjgHwMhpcSU5lNURyEioteMCxuHozeOYkW7FaqjkIkz+Iy3EGIBgEMA/gQQLqVMffZ4fQBXpZQRGQmTGWy8SZ8iYyNRaXYlnOl9Bh9/9LHqOERE9JqEpASU+KUENnpvRKWClVTHIROmYsb7CoBoAL4AQoUQIUKIAABWANpmcJ9kprQw5zZq9yj0qtKLTfdbaKF+9G6sn3axdq/KmSUnRtQdgRE7RqiOkiasn3myyshGUspRzz5dAwBCiOwAqgOoCeCybqIRGYdzd85hw4UNuNjvouooRET0Hl9W/hIT90/Ensg9+MTxE9VxiN6g0xlvIYSLlFJJ481RE9KXdr+3Q5WCVTC0zlDVUYiI6AMWhi/E7COzsbfrXr7JGemF0tsJCiF+FkIMFkJUA5AihOie2X0SGYujN47iz6g/0a96P9VRiIgoDXxcfRCbGIvNlzarjkL0Bl3cTvAsgIsAWgD4FYCrDvZJZsSY59y+3/k9htcdjhzWOVRHMVrGXD/6MNZPu1i7t7O0sMQoj1H4fuf3MOa/hLN+5kkXjbeQUq6XUo6QUtYHsEcH+yRS7uD1gzjx9wl0r8Q/4hARaUmrMq3wJPUJ1p9frzoK0SsyPeMthGgBwA/ARgCnATSUUo7VQbb05uCMN+lU08VN0bp0a/So0kN1FCIiSqf159dj5M6RONbjGCwE3y+QdEfpjLeUch2AbwE4AegBYEtm90mkWlhUGC7cvYCuFbuqjkJERBngVdILWSyzYNWZVaqjEL2QocZbCLHppc8HAmgNYC6A/wEoqJtoZC6Mcc7t+53fY+QnI5HFMovqKEbPGOtHacf6aRdr935CCIz2GI3A3YFISU1RHecNrJ95ylDjLaX89KXFeADHAYwG8COeXmSZJkKIZkKIc0KIC0KIN+7VJoSwEUKsF0IcF0KcFEL4ZyQvUXrsjNiJqLgo+Ln5qY5CRESZ0Kx4M9hmtcVvp39THYUIgG5mvJ0BOEgpw4QQWQBYSSkfpWE7CwAXADQEEIOnb0HfUUp57qV1vgNgI6X8TgiRF8B5AB9LKZ+8ZX+c8aZMk1Ki7q910aNyD3R266w6DhERZdK2K9vQO6Q3zvQ5AyuLDL1vINErlM54A6gM4Nqzz3Onpel+phqAi1LKSCllMoDlePNsuQSQ69nnuQDcfVvTTaQrW69sxd3Hd+Ht6q06ChER6UBD54YomKsglpxYojoKkU4abxf82xzbPLvLSVo44N+GHQCuP3vsZf8DUFYIEQMgHMCAzAQl42Qsc25SSozcORKB9QJhaWGpOo5mGEv9KGNYP+1i7dLm+az36D2jkZySrDrOC6yfedJF450EwFcI0RzAHQDZdbDP55oCOCalLASgIoBpQoiPdLh/ohdCLobgUfIjtCvXTnUUIiLSoXpO9eBs54zg48Gqo5CZ08WwUyqANQBqA+gD4HAat4sGUPSl5cLPHntZVwBjAUBKeVkIEQGg9Lu+hr+/P5ycnAAAdnZ2cHd3h4eHB4B/f7PksvEte3h4KM+zc+dOfL3xa4ztNhYWwkJ5Hi0tG0P9uJzxZdaPy+ayHFQ/CB1XdYRTrBOsLa2V5+GydpaPHz+O2NhYAMDVq1eRGbq4uLITgA1SynghREEAn0spZ6VhO0s8vViyIYAbAP4C0ElKefaldaYBuCWlHCWE+BhPG243KeW9t+yPF1dShm28sBEjdozgGy0QEZmwpoubom2Ztviy8peqo5CGqb64cjmA4s8+LwAgT1o2klKmAOgLIBRP3/FyuZTyrBCihxDiq2er/QCglhDiBICtAL59W9NN2vb8t0tVpJQYtXsURn4ykk13BqiuH2UO66ddrF36ff/J9/gx7EejmPVm/cxTpkdNnp1mPv5s8YGU8sd0bPsHgFKvPTbrpc9v4OmcN5HebL60GYlPEtG6TGvVUYiISI9qF60NF3sXLAxfiG6VuqmOQ2ZIF6MmP+PpHUn2ArgFoJGUcq4OsqU3B0dNKN2klKgxrwYG1xyM9uXaq45DRER6tjdyL7qs7YLzfc/D2tJadRzSIL2OmgghGn5glbMALuLpPbh/BeCakSBEKmy5vAXxSfFoW7at6ihERGQAdR3rwsnOCYtPLFYdhcxQWgZaR33geSGlXC+lHCGlrA9gjw5ykRkx1JxbREQkfH1HoX79APj6jsKVK1c5260DnFPUNtZPu1i7jAuoF4Axe8fgSeqTN342REREGiQD62ee0jLjXU0I8QuA/0gp497y/E0hxCoAG/H0IsmSugxIpAsREZFo3PgXXL48CkBOAAnYEdkFOdrcRruyvG83EZE5qedUDw42DpiyYypm9Ix55WfDgQMB2Lq1H5ydHVXHJBP0wRlvIURvAIsA9ANwXUq58C3ruADwA2APIFhKeVQPWd+LM970Pr6+o7BkyRA8/Y8VACTwRW3UylIQf85YpTIaEREpsCNiB1rObYeH464AqbYvPZMAH5+JWLw4QFk2Mm56nfGWUk6XUj58dreSY0KICUKIiq+tc1lKGSCl7K+i6Sb6kOjoVPzbdANw3gHkuIss58sqy0REROrUd6oPi0c5gfIbXnsmJ2JiUpVkItOXrsFWKeVJKeU3AMoJIQKEEHZ6ykVmxBBzbg4OFgAS/n2g3mhgzzdwKKSLN281b5xT1DbWT7tYu8wRQqByvAfwyWhApLz0TAIKFdL/dT+sn3nK0HeWlHIxgMkABgohuug2EpHuBQX5w8UlAEAC4LQLyBWNYo/OICjIX2kuIiJSZ853o5E1NR4o93yKNgEuLgH82UB6k5YZ7+GvvymOEEIAcADgCKA1gJoA+kgpj+kr6Idwxps+JCIiEiNHBmNT3gUolVQBS7+ZwotniIjMXHDYAgzc9A0qHugJh0KWCAry588Geq/MzHinpfE+j6dvC+/40ocD/r0jSjKevoHOaSlli4yE0AU23pQWB68fRPuV7XGp3yW+cQIREUFKiapzquI/n/wHLUu3VB2HNECvF1cCKAagJYDceHq7wBkAOgOoDaAwgOxSyuIqm27SNkPOuY3ZOwbf1vqWTbcOcU5R21g/7WLtdEMIgRF1R2DM3jEw5Ak81s88paXxXiqldJNSfi6l7CulnCClXCGlPCClvMHTzKQV4TfDcTjmML6o+IXqKEREZERalG6Bx8mPEXo5VHUUMnFpGTXJJaV8aKA8GcZRE/qQDis7oGqhqhhSa4jqKEREZGSWnFiCWUdmYU9XvgE3vZ++7+Nt9E030Yecv3MeOyN2omeVnqqjEBGREepQvgOiH0Zjb+Re1VHIhOn/RpVEH2CIObdxf45Dv2r98FGWj/T+tcwN5xS1jfXTLtZOt6wsrPBdne8wZu8Yg3w91s88sfEmk3c19irWn1+PvtX6qo5CRERGzM/ND6dvn8bhmMOqo5CJ+uCMt1ZwxpvepXdIb9hmtcXYRmNVRyEiIiM39eBU7Ly6E2s6rFEdhYyUXu/jrRVsvOltbjy8gXLTy+Fc33PInzO/6jhERGTkHiU/QrEpxbDNbxvK5y+vOg4ZIX3fx5tIr/Q55zZp/yT4ufmx6dYjzilqG+unXaydfuSwzoFBNQZhbJh+/0rK+pknNt5ksu4+uotfj//K2wcSEVG69KraC6GXQ3H53mXVUcjEcNSETNbo3aMRFReFuZ/PVR2FiIg0ZuSOkbjz6A5meM5QHYWMDGe8wcabXpWQlADnKc7Y23UvSuUtpToOERFpzO2E2yj1v1I40+cMCnxUQHUcMiKc8SZN08ec27xj81DXsS6bbgPgnKK2sX7axdrpV76c+eDt6o0pB6boZf+sn3li400mJzklGZP2T8LQ2kNVRyEiIg0bXHMw5hydg7jEONVRyERw1IRMzqLwRQgOD8Z2v+2qoxARkcb5rvZFhY8r4Nva36qOQkaCoyZEz6TKVPz3z/9iWO1hqqMQEZEJGFp7KCYfmIzEJ4mqo5AJYONNyulyzi3kQgiyWGZBo2KNdLZPej/OKWob66ddrJ1huH7siooFK2Jh+EKd7pf1M09svMmk/PfP/2JYnWEQIkN/ASIiInrDsNrDMGHfBKSkpqiOQhrHGW8yGWFRYfBf64/zfc/D0sJSdRwiIjIRUkrU+bUOBlYfiHbl2qmOQ4pxxpsIwLiwcfim1jdsuomISKeEEBhaeyjG/TkOPMlHmcHGm5TTxZzbyb9P4uiNo+ji3iXzgShdOKeobayfdrF2huVZ0hOJTxKxPUI3d8xi/cwTG28yCeP3jUf/6v2RzSqb6ihERGSCLLQAfxkAACAASURBVITF07PeYeNURyEN44w3ad61uGtwn+WOy/0vwy6bneo4RERkopJSkuAy1QXrOq5DpYKVVMchRTjjTWZtysEp8HfzZ9NNRER6lcUyCwZUH4BJ+yepjkIaxcablMvMnFtcYhx+Pf4rBtQYoLtAlC6cU9Q21k+7WDs1vqz0Jf649AciYyMztR/Wzzyx8SZNm31kNpoXb46itkVVRyEiIjNgm80WX7h/gSkHp6iOQhrEGW/SrKSUJBSbUgwbvTfCvYC76jhERGQmrsVdg9tMN1wZcIVjjmaIM95klpadXIYy+cqw6SYiIoMqYlsEniU9MfPwTNVRSGPYeJNyGZlzk1Ji4v6JGFJziO4DUbpwTlHbWD/tYu3UGlxzMKYenIp/nvyToe1ZP/PExps0acvlLRAQaOLSRHUUIiIyQ24F3OD6sSuWnlyqOgppCGe8SZMaLWwEPzc/+Ln5qY5CRERmauvlrRi4ZSBO9ToFITI08ksaxBlvMitHbxzFuTvn0LF8R9VRiIjIjDUq1gjWFtbYfGmz6iikEWy8Sbn0zrlN2j8JA6oPQBbLLPoJROnCOUVtY/20i7VTTwiBIbWGYOK+ienelvUzT2y8SVMiYyPxx6U/8FXlr1RHISIiQodyHXDp3iUcjjmsOgppAGe8SVMGbxkMAJjUlG/XS0RExmHSvkk4cuMIlrbhhZbmIDMz3my8STMe/PMAzlOccfSro3C0c1Qdh4iICAAQlxgH5ynOCO8ZjiK2RVTHIT3jxZWkaWmdc5t/bD4aFWvEptvIcE5R21g/7WLtjIdtNlv4ufnhf3/9L83bsH7mSWnjLYRoJoQ4J4S4IIQY+o51PIQQx4QQp4QQOw2dkYxDSmoKph6cikE1BqmOQkRE9Ib+1ftj3rF5iE+KVx2FjJiyURMhhAWACwAaAogBcAhARynluZfWsQWwD0ATKWW0ECKvlPLOO/bHURMTtvrsakzcNxH7uu1THYWIiOit2qxogwZODdCnWh/VUUiPtDpqUg3ARSllpJQyGcByAC1eW8cbwCopZTQAvKvpJtP30/6feLabiIiM2qAagzD54GSkylTVUchIqWy8HQBce2n5+rPHXlYSQG4hxE4hxCEhRGeDpSOD+dCc21/Rf+H6g+toVaaVYQJRunBOUdtYv/+3d+fhVVbn+sfvJ4ACigJFKyjzqAwGBJSqiLQo8rNqKSCTEFFBFLHa9mjbc7Stv/bYczwOoCIcQEDA4AyOSAVnRiHMyBhkEBAxVBAkJOv8sTc2pCEDSfbaa+/v57q4yNqsvLnDk+HJm+d9d7ioXfy5tO6lqlG5ht5Y/0aRe6lfcqroO0ARKkpqJ6mrpNMkzTez+c65jQVtTktLU4MGDSRJ1atXV2pqqrp06SLpnx/grMNbP7bgMfWo1EMff/hxXORhzZo163hYHxMveVh3kZnp6gpX64FnH9B1D19X6P5j4ik/64LXGRkZysrKkiRlZmaqNHzOeF8i6Y/Oue7R9f2SnHPub3n23CepsnPuT9H1eElvO+deLuB4zHgnoG37tyl1bKo2j9ysMyuf6TsOAACFys7JVqNRjTSr7yy1rd3WdxyUg1BnvBdLamJm9c3sFEl9Jc3Kt2empMvMrIKZVZV0saS1Mc4Jj55c9KQGtRlE0w0ACEKlCpV0V8e79NiCx3xHQRzy1ng753IkjZD0rqTVktKdc2vNbJiZDY3uWSdptqQVkhZIGuecW+MrM8pH/l+7HXPgyAFNWDZBIy8eGdtAKJET1Q9hoH7honbx67Z2t+mN9W9o57c7T7iH+iUnrzPezrl3JDXP99jYfOtHJD0Sy1yID5MyJqlLgy5qWKOh7ygAABRbjSo1NKD1AD216Cn95ad/8R0HcYSnjEdcysnNUfMnm2vyDZN1ab1LfccBAKBENu7bqE4TOmnrr7aqaqWqvuOgDIU64w2c0Fsb3lKNKjX0k7o/8R0FAIASa1KziX5S9yeaumKq7yiIIzTe8K6gObdRi0bp7ovvltlJ/UCJGGJOMWzUL1zULv6N7DhSoxaOUkG/kad+yYnGG3Fn9Z7VWrVnlXpf0Nt3FAAATlrXhl0lSXO3zPWcBPGCGW/EndvfuF21T6+tB7s86DsKAAClMu6zcXpzw5ua2Xem7ygoI8x4I2HsO7RPM1bP0LD2w3xHAQCg1Aa2GahPt32qTfs2+Y6COEDjDe/yzrlNWDpBP2/2c51z+jn+AqFEmFMMG/ULF7ULQ9VKVTUkdYieWvzUcY9Tv+RE4424cTT3qJ5c/CRPmAMASCh3dLhDk5dP1rfff+s7Cjxjxhtx4+U1L+vRBY/qkyGf+I4CAECZ6vVCL3Vp0EUjOo7wHQWlxIw3EsKoRaM0siNnuwEAiWfkxSM1etFo5bpc31HgEY03vHv//feVsStDm/ZtUs/ze/qOgxJiTjFs1C9c1C4sl9e7XFUrVdXsjbMlUb9kReONuDBq4Sjd0eEOVapQyXcUAADKnJlFnlBn0SjfUeARM97w7quDX6nZk8204a4NqlW1lu84AACUi8NHD6v+4/X1QdoHalGrhe84OEnMeCNo4z4bp54tetJ0AwASWuWKlXVbu9s0euFo31HgCY03vDqae1SPpz+uuy6+y3cUnCTmFMNG/cJF7cI0vP1wTV81XW+++6bvKPCAxhtezVw3U+dUO0ep56T6jgIAQLk794xz1a1RN7276V3fUeABM97w6srJV2rYRcPUt1Vf31EAAIiJDzI/0LA3hmntnWtldlKjwvCIGW8EadWeVVq3dx23EAQAJJXO9TurUoVKem/Le76jIMZovOHN04uf1rCLhunTjz71HQWlwJxp2KhfuKhduMxM3VK66anFT/mOghij8YYX+w/v1/OrntfQi4b6jgIAQMx1a9RNH279UFuztvqOghhixhtejFo4Sp9u+1TpvdJ9RwEAwIt73rlHlStW1n/+7D99R0EJMOONoOS6XD21+CmN6DjCdxQAALy5o8MdmrBsgg4fPew7CmKExhsx997m91S5YmVdWvdSScwpho76hY36hYvahe39999X0x81Vbva7fTC6hd8x0GM0Hgj5p5c/KRGdBjBLZQAAElvRMcRenLRk75jIEaY8UZMZWZlqv249tr6q6067ZTTfMcBAMCrnNwcNR3dVOm90tXx3I6+46AYmPFGMMYsHqPBFw6m6QYAQFKFlAq6o8MdnPVOEjTeiJlD2Yf0bMazGt5h+HGPM6cYNuoXNuoXLmoXtrz1G9J2iF5f/7q+OviVv0CICRpvxMwLq19Q+zrt1aRmE99RAACIGzWr1FTPFj01YdkE31FQzpjxRsxcMv4S/eHyP+jnzX/uOwoAAHFlyc4l6v1ib228a6MqpFTwHQeFYMYbcW/pl0v15YEv1aNpD99RAACIO+3rtFetqrU0e9Ns31FQjmi8ERNjFo/R0HZDC/wpnjnFsFG/sFG/cFG7sBVUv+Hth+vpxU/HPgxihsYb5W7/4f16ae1LuqXdLb6jAAAQt/q26qv52+crMyvTdxSUE2a8Ue5GLxytT7Z9ovRe6b6jAAAQ1+555x5VqVRFf/3pX31HwQkw44245ZzTmCVjNLz98KI3AwCQ5G5vf7smLpuoIzlHfEdBOaDxRrn6cOuHMjN1rt/5hHuYUwwb9Qsb9QsXtQvbierXvFZztTy7pV5Z+0psAyEmaLxRrsYsGaPbL7pdZif1GxkAAJLO8PbDNWbJGN8xUA6Y8Ua52X1gt1o81UKZd2fqzMpn+o4DAEAQsnOyVf/x+ppz0xy1PLul7zjIhxlvxKUJyyao1/m9aLoBACiBShUq6dZ2t+qZJc/4joIyRuONcpGTm6Oxn43V8A5FX1TJnGLYqF/YqF+4qF3Yiqrfbe1u07SV03TgyIHYBEJM0HijXLy14S2dc/o5ale7ne8oAAAEp+6ZddW5fmdNXznddxSUIWa8US56TOuhPi37KC01zXcUAACCNHvjbN3/3v1aOnQpNymII8x4I65kZmVq0Y5FurHljb6jAAAQrG6Nu+nb77/V4p2LfUdBGaHxRpkbv3S8BrYZqCqVqhRrP3OKYaN+YaN+4aJ2YStO/VIsRbe1u01jl4wt/0CICRpvlKnsnGxNXDZRQy8a6jsKAADBS0tN0yvrXtH+w/t9R0EZYMYbZerVta/q0QWP6qObP/IdBQCAhNDnxT7q0qCL7uhwh+8oEDPeiCNjPxurYRcN8x0DAICEMeyiYRr72VhxgjF8NN4oM1u+2aIlO5fol+f/skSvx5xi2Khf2KhfuKhd2EpSvysbXqmDRw5q0Y5F5RcIMeG18Taz7ma2zszWm9l9hezrYGbZZtYzlvlQMuOXjtdNbW4q9kWVAACgaCmWoqEXDdXYz7jIMnTeZrzNLEXSekk/lbRT0mJJfZ1z6wrYN0fSIUkTnXOvnOB4zHh7lJ2TrXqP19PcQXN1/lnn+44DAEBC2XNwj5o/2Vxb7t6i6pWr+46T1EKd8e4oaYNzbqtzLltSuqTrC9h3l6SXJO2JZTiUzKzPZ6lpzaY03QAAlIOzTztbVzW+SlNXTPUdBaXgs/E+V9K2POvt0cd+YGZ1JN3gnBsjiadsimPjlo476YsqmVMMG/ULG/ULF7UL28nUb9hFwzTus3FcZBmwir4DFOFxSXlnvwttvtPS0tSgQQNJUvXq1ZWamqouXbpI+ucHOOuyX2/+ZrMWfLxAv679ax0TT/lYs2bNOhHXx8RLHtYlWx9Tkte/ssGV+nrt13r6xad1Z5874+r9SeR1RkaGsrKyJEmZmZkqDZ8z3pdI+qNzrnt0fb8k55z7W549m4+9KKmWpIOShjrnZhVwPGa8Pfnd33+n73O+16NXP+o7CgAACe2/P/lvrf5qtSbdMMl3lKRVmhlvn413BUmfK3Jx5ZeSFknq55xbe4L9z0p6nYsr48uRnCOq91g9vZ/2vlrUauE7DgAACe2rg1+p6eim2nL3FtWoUsN3nKQU5MWVzrkcSSMkvStptaR059xaMxtmZkMLepWYBkSxzPp8llrUalGqpjv/r90QFuoXNuoXLmoXtpOt31mnnaVrml7DRZaB8jrj7Zx7R1LzfI8VeJNK59yQmIRCifzv0v/Vbe1u8x0DAICkMbTdUI18Z6RGdBwhM+49ERJvoyZljVGT2MvMylT7ce21/d7tqlyxsu84AAAkhVyXq2ajm2n6L6er47kdfcdJOkGOmiB8zy57Vv1b96fpBgAghlIsRbe0vUXjl473HQUlROONk5KTm6OJGRN1a7tbS30s5hTDRv3CRv3CRe3CVtr6paWm6cU1L+rAkQNlEwgxQeONkzJ702zVqVZHbX7cxncUAACSTu1qtXVF/Ss0Y9UM31FQAsx446T0nNFT1zS5RrddxIWVAAD48Mb6N/SXj/6i+bfM9x0lqTDjjZjadWCX5mXOU99WfX1HAQAgaXVv0l1f7P9Cq/as8h0FxUTjjRKbnDFZPVv0VLVTq5XJ8ZhTDBv1Cxv1Cxe1C1tZ1K9iSkXdnHozF1kGhMYbJeKc0/hl48vkokoAAFA6Q9oO0dQVU3X46GHfUVAMzHijRD7I/EB3vnWnVg5fyU37AQCIA92e66Zb2t7CCGiMMOONmDl2tpumGwCA+HBr21sZNwkEjTeK7ZtD3+j1z1/XwDYDy/S4zCmGjfqFjfqFi9qFrSzrd0OLG7R893Jt2repzI6J8kHjjWKbtnKaujfprlpVa/mOAgAAok6teKoGth6oicsm+o6CIjDjjWJxzil1bKr+56r/0c8a/cx3HAAAkMfqPavV7blu+uKeL1QxpaLvOAmNGW+Uu8++/Ez/+P4f6tqwq+8oAAAgn5Znt1T96vX19oa3fUdBIWi8USwTl03UkNQhSrGy/5BhTjFs1C9s1C9c1C5s5VG/W9reookZjJvEMxpvFOlQ9iHNWD1Dg1MH+44CAABOoE/LPpq3ZZ72HNzjOwpOgBlvFOn5lc9r0vJJmj1wtu8oAACgEGmvpanNj9vo3k73+o6SsJjxRrmamBEZMwEAAPFtSNshmrBsgjgZGZ9ovFGozKxMLftyma5vcX25vQ3mFMNG/cJG/cJF7cJWXvW7vN7lOpJzRIt3Li6X46N0aLxRqMkZk9WvVT9VrljZdxQAAFAEM9PNqTdzT+84xYw3TijX5arRE4306o2vqm3ttr7jAACAYtj+j+1qM6aNtt+7XVUrVfUdJ+Ew441yMW/LPNWoUoOmGwCAgJx3xnm65LxL9MraV3xHQT403vjBli1bNXDgn3TllQ9q4MA/adTHo2NyUSVzimGjfmGjfuGidmEr7/oNaTtEE5dN/Jfv7Vu2bC3Xt4vC8ZyikBRpurt1G61Nm/4k6TSp8nal1G+iBzv80Xc0AABQQj9v9nMNmzVMXX7xkL5Y/oSk0yQd1IIFD2rOnLvUsGF93xGTEjPekCQNHPgnTZv2G0U+MSW1HyM1+LsGnNpGU6c+6DUbAAAouRa/ukSfr+gizXs4z6MHNWDAI3xvLwVmvFFqO3bk6oemW5LaTpSWDdXOnbneMgEAgJNXbeOFUup0yXLyPHoa39s9ovGGJOncc1MkHYwsfrxCOn2XtLmT6tQp/w8R5hTDRv3CRv3CRe3CFov6Na9eRzpYS2r0Xp5HD8bkezsKxv88JEkPPZSmxo0flHRQSn1Wyuivxo3+rIceSvMbDAAAnJSHHkpTrW01pbZjo48cVOPGD/K93SNmvPGDLVu26vf/MV4v139U1+y8VY8/cC8XXwAAELCMdSvUcXpHXbxgpOqffZoeeiiN7+2lVJoZbxpvHOe1da/p0fmP6sObP/QdBQAAlIEbX7pRXep30fAOw31HSQhcXIkyM3n5ZKWlpsX0bTKnGDbqFzbqFy5qF7ZY1i/twjRNXj45Zm8PJ0bjjR98dfArzdsyT70v6O07CgAAKCPdGnfTF/u/0Nqv1vqOkvQYNcEPnljwhD778jNN+cUU31EAAEAZum/OfTIzPfyzh4vejEIxaoIyMWn5JA2+cLDvGAAAoIwNTh2s51Y8p5zcnKI3o9zQeEOStHzXcn393de6suGVMX/bzCmGjfqFjfqFi9qFLdb1u+CsC3RutXM1Z/OcmL5dHI/GG5IiF1UOvnCwUowPCQAAElFaKhdZ+saMN5Sdk63zHjtPnwz5RE1qNvEdBwAAlIN9h/ap0RONlPmrTFWvXN13nGAx441SeXvj22pasylNNwAACaxmlZrq1ribZqya4TtK0qLxhiZlTIr5vbvzYk4xbNQvbNQvXNQubL7ql3ZhmiYtn+TlbYPGO+nt/W6v5m6Zy727AQBIAlc3uVqZWZn6fO/nvqMkJWa8k9zohaO1cMdCTe051XcUAAAQA79997eqVKGS/vrTv/qOEiRmvHHSuHc3AADJZXDqYE1ZPoV7entA453EVuxeoT0H96hrw65eczCnGDbqFzbqFy5qFzaf9Wt1diudc/o5em/Le94yJCsa7yQ2ZfkU3dTmJlVIqeA7CgAAiCHu6e0HM95J6mjuUdV7rJ7mDp6rFrVa+I4DAABiaO93e9VkVBNtu2ebqp1azXecoDDjjRJ7b/N7Ou+M82i6AQBIQrWq1tIVDa7Qy2tf9h0lqdB4J6kpK6Zo0IWDfMeQxJxi6Khf2KhfuKhd2OKhfoPaDNKU5VN8x0gqXhtvM+tuZuvMbL2Z3VfAv/c3s+XRPx+bWWsfORPNt99/qzfXv6m+rfr6jgIAADy5ttm1Wr57ub7Y/4XvKEnD24y3maVIWi/pp5J2Slosqa9zbl2ePZdIWuuc229m3SX90Tl3yQmOx4x3MU3KmKRX172qmX1n+o4CAAA8uv2N21X/zPr63eW/8x0lGKHOeHeUtME5t9U5ly0pXdL1eTc45xY45/ZHlwsknRvjjAlpyvIpGtQmPsZMAACAP4MuHKQpK6aIk5ex4bPxPlfStjzr7Sq8sb5V0tvlmigJbM3aquW7l+vaZtf6jvKDeJhzw8mjfmGjfuGidmGLl/p1Oq+TsnOytWTnEt9RkkJF3wGKw8yulHSzpMsK25eWlqYGDRpIkqpXr67U1FR16dJF0j8/wJN9/WmFT9Xngj6a//H8uMjDmjVr1qxPbn1MvORhXbL1Mb7zfPDBB7o893JNWT5FHc7t4D1PPK4zMjKUlZUlScrMzFRp+JzxvkSRme3u0fX9kpxz7m/59rWR9LKk7s65TYUcjxnvIjjndP5T52vi9RP1k7o/8R0HAADEgc3fbNbF4y/Wjnt36JQKp/iOE/dCnfFeLKmJmdU3s1Mk9ZU0K+8GM6unSNN9U2FNN4pnyc4lOpp7VJ3O6+Q7CgAAiBONajRSi1ot9M7Gd3xHSXjeGm/nXI6kEZLelbRaUrpzbq2ZDTOzodFt/yGppqSnzWyZmS3yFDchTFkeuXe32Un9kFZu8v/aDWGhfmGjfuGidmGLt/pxT+/Y8Drj7Zx7R1LzfI+NzfPybZJui3WuRHQk54jSV6dr4a0LfUcBAABxpnfL3vrNnN9o36F9qlmlpu84CcvbjHdZY8a7cDPXzdQj8x/RRzd/5DsKAACIQ31e7KOuDbvq9va3+44S10Kd8UYMPbfiOe7dDQAATmjQhYP03IrnfMdIaDTeSeCbQ99ozuY56t2yt+8oBYq3OTeUDPULG/ULF7ULWzzW7+rGV2vjvo3auG+j7ygJi8Y7Cby45kVd1fgqVa9c3XcUAAAQpypVqKQbW96oaSum+Y6SsJjxTgKdn+2sX3f6ta5vcb3vKAAAII4t2rFIA14ZoPUj1sfdXdDiBTPeOKGtWVu15qs1uqbpNb6jAACAONehTgeZTIt3LvYdJSHReCe46Sunq9cFveL6majicc4NxUf9wkb9wkXtwhav9TMzDWg9QFNXTPUdJSHReCcw55ymrpyqgW0G+o4CAAACMaDNAM1YPUPZOdm+oyQcZrwTWMauDN2QfoM2371ZKcbPWAAAoHg6TeikBzo/wKhqAZjxRoGmrpiqAa0H0HQDAIASGdh6oKauZNykrNGRJaic3BxNXzldA9oM8B2lSPE654bioX5ho37honZhi/f69WnZR2+uf1MHjhzwHSWh0HgnqPcz31ftarV1wVkX+I4CAAACc9ZpZ+myepfptXWv+Y6SUJjxTlA3z7xZrc9urXs73es7CgAACFD6qnRNypikdwa+4ztKXGHGG8f5Lvs7vbbuNfVt1dd3FAAAEKjrml+nhTsWateBXb6jJAwa7wT0+uevq0OdDqpTrY7vKMUS73NuKBz1Cxv1Cxe1C1sI9ataqaqua36d0lel+46SMGi8E9C0ldM0oHX8X1QJAADi28DWAzVt5TTfMRIGM94JZu93e9V4VGNtu2ebzjj1DN9xAABAwHJyc1T3sbqaO3iuWtRq4TtOXGDGGz94YfUL6tG0B003AAAotQopFdSvVT9NW8FZ77JA451gQhwzCWHODSdG/cJG/cJF7cIWUv0GtBmgaSunicmC0qPxTiCZWZn6fO/nuqrxVb6jAACABNH2nLY6pcIpWrhjoe8owWPGO4E8/PHDyszK1DPXPuM7CgAASCB//uDP+vq7r/XENU/4juIdM96QJE1fOV39W/f3HQMAACSYfq36acbqGTqae9R3lKDReCeIlbtX6pvD3+iyepf5jlJiIc254V9Rv7BRv3BRu7CFVr+mP2qqumfW1bwt83xHCRqNd4J4ftXz6teqn1KMkgIAgLLXv1V/TV813XeMoDHjnQCcc2o0qpFe6fOK2tZu6zsOAABIQDu/3amWT7fUl7/+UpUrVvYdxxtmvJPcgu0LVLliZaWek+o7CgAASFB1qtVRu9rt9NaGt3xHCRaNdwKYvnK6+rfqL7OT+uHLu9Dm3HA86hc26hcuahe2UOvXv1V/TV/JuMnJovEO3NHco3phzQvq17qf7ygAACDB9Ty/p+ZsnqP9h/f7jhIkZrwDN3vjbD3w/gNaeCs3tQcAAOXvhvQbdEOLG5SWmuY7ihfMeCex51c9r/6tuHc3AACIjf6t++v5Vc/7jhEkGu+AHco+pJmfz1Sfln18RymVUOfcEEH9wkb9wkXtwhZy/a5tdq0W7Vik3Qd2+44SHBrvgL254U1dVPsi1a5W23cUAACQJKpWqqprm12rF1a/4DtKcJjxDljPGT11bbNrNaTtEN9RAABAEnl7w9v684d/1vxb5vuOEnPMeCehrMNZem/Le+p5fk/fUQAAQJL5WaOfadO+Tdr8zWbfUYJC4x2oV9e+qq4Nu6p65eq+o5RayHNuoH6ho37honZhC71+lSpUUq8Leun5lVxkWRI03oFKX52ufq24dzcAAPCjX6t+mrF6hu8YQWHGO0B7Du5Rs9HNtOPeHTrtlNN8xwEAAEko1+Wq3mP1NHvgbLU8u6XvODHDjHeSeXnNy+rRtAdNNwAA8CbFUnRjyxs5610CNN4BSl+drr6t+vqOUWZCn3NLdtQvbNQvXNQubIlSv76t+ip9VbqSZeqgtGi8A7P9H9u1cvdKXd34at9RAABAkmtfp71yXI6W7VrmO0oQmPEOzGPzH9PKPSs18fqJvqMAAADoD+/9Qdm52fqvbv/lO0pMMOOdRBJtzAQAAIStb6u+mrF6hnJdru8ocY/GOyCbv9msLd9sUdeGXX1HKVOJMueWrKhf2KhfuKhd2BKpfq3ObqXTTzldC7Yv8B0l7tF4B2TGqhnqdUEvVUyp6DsKAACApMjoRd+WkYssUThmvANy4TMXavQ1o9W5fmffUQAAAH6w/uv1umLSFdp+z3ZVSKngO065YsY7Caz5ao32frdXl9W7zHcUAACA4zT7UTPVqVZHH2z9wHeUuEbjHYgZq2boxpY3KsUSr2SJNOeWjKhf2KhfuKhd2BKxfoybFM1rF2dm3c1snZmtN7P7TrBnlJltMLMMM0uNdcZ44JxL6LuZZGRk+I6AUqB+YaN+4aJ2YUvE+vVp2UevrH1FR3KO1J49iQAACSVJREFU+I4St7w13maWIulJSVdLaimpn5m1yLfnGkmNnXNNJQ2T9EzMg8aBjF0Zys7JVoc6HXxHKRdZWVm+I6AUqF/YqF+4qF3YErF+9avXV/NazfX3zX/3HSVu+Tzj3VHSBufcVudctqR0Sdfn23O9pCmS5JxbKOlMM/txbGP6l74qcrbb7KTm+AEAAGKCcZPC+Wy8z5W0Lc96e/SxwvbsKGBPQnPOacbqyHx3osrMzPQdAaVA/cJG/cJF7cKWqPXr3bK3Xl//ug5lH/IdJS55u52gmf1S0tXOuaHR9UBJHZ1zI/PseV3SfzrnPo2u/y7p35xzSws4XmLfSxAAAABx4WRvJ+jzmVh2SKqXZ31e9LH8e+oWsUfSyf8HAAAAALHgc9RksaQmZlbfzE6R1FfSrHx7ZkkaJElmdomkLOfc7tjGBAAAAErP2xlv51yOmY2Q9K4iPwBMcM6tNbNhkX9245xzb5lZDzPbKOmgpJt95QUAAABKI2GeMh4AAACIZ8E+DaKZ1TCzd83sczObbWZnnmDfmWb2opmtNbPVZnZxrLPieMWtXXRvipktNbP8Y0jwpDj1M7PzzGxu9HNupZmNLOhYiA2erCxsRdXPzPqb2fLon4/NrLWPnPhXxfnci+7rYGbZZtYzlvlQuGJ+7exiZsvMbJWZzSvqmME23pLul/R351xzSXMl/e4E+56Q9JZz7nxJF0paG6N8OLHi1k6S7pa0JiapUFzFqd9RSfc651pK6iTpzvxPkIXY4MnKwlac+knaLKmzc+5CSf9f0v/GNiUKUszaHdv3sKTZsU2IwhTza+eZkp6SdK1zrpWk3kUdN+TG+3pJk6MvT5Z0Q/4NZnaGpMudc89KknPuqHPuH7GLiBMosnZS5KyppB6SxscoF4qnyPo553Y55zKiLx9Q5AfepLoHfxzhycrCVmT9nHMLnHP7o8sF4nMtXhTnc0+S7pL0kqQ9sQyHIhWnfv0lveyc2yFJzrm9RR005Mb77GN3OHHO7ZJ0dgF7Gkraa2bPRscVxplZlZimREGKUztJekzSbyVxIUJ8KW79JElm1kBSqqSF5Z4MBeHJysJWnPrldaukt8s1EYqryNqZWR1JNzjnxkjitsjxpTife80k1TSzeWa22MxuKuqgPu/jXSQzmyMp71kXU6QJ+/cCthfUnFWU1E7Snc65JWb2uCK/Jn+wrLPieKWtnZn9P0m7nXMZZtZFfEGKqTL43Dt2nNMVOZNzd/TMN4ByYmZXKnL3r8t8Z0GxPS4p7+ww3+vCcqzP7CrpNEnzzWy+c25jYa8Qt5xz3U70b2a228x+7JzbbWbnqOBf0WyXtM05tyS6fknHf4CjnJRB7S6VdJ2Z9ZBURVI1M5vinBtUTpGRRxnUT2ZWUZHPueecczPLKSqKVqZPVoaYK079ZGZtJI2T1N05902MsqFwxalde0npZmaSakm6xsyynXPcUMC/4tRvu6S9zrnDkg6b2YeKXE94wsY75FGTWZLSoi8PlvQv39ijvw7fZmbNog/9VFyoFw+KU7vfO+fqOecaKfLkSnNpuuNGkfWLmihpjXPuiViEwgnxZGVhK7J+ZlZP0suSbnLObfKQEQUrsnbOuUbRPw0VOVFxB0133CjO186Zki4zswpmVlXSxSriJh4hN95/k9TNzD5XpKF+WJLMrLaZvZFn30hJ08wsQ5GfQv4a86TIr7i1Q3wqsn5mdqmkAZK6Rm+ztNTMuntLnMScczmSjj1Z2WpJ6ceerMzMhkb3vCVpi0WerGyspDu8BcZxilM/Sf8hqaakp6Ofb4s8xUUexazdca8S04AoVDG/dq5T5G40KxS5sHmcc67QE7w8gQ4AAAAQAyGf8QYAAACCQeMNAAAAxACNNwAAABADNN4AAABADNB4AwAAADFA4w0AAADEAI03AMQxM3vQzHLz/PnSzF43s9ZldPxHzGxLnvVgM8uJPhkEAKAM0XgDQPzLUuQZ0S6RdLekZpLeNbPqZXBsp+OfuOMNSZ2cc9+VwbEBAHlU9B0AAFCko865xdGXF5nZVknzJXWXlF6Wb8g597Wkr8vymACACM54A0B4lkf/ritJZlbVzEab2TozO2hmm83sSTOrlveVzOxMM5tuZt+a2Q4z+33+A5tZWnSkpWp0fUV0fUG+ffPM7IU86wvM7G0z+9rMDpjZGjMbXubvOQAEjDPeABCe+tG/N0f/rqrI1/N/l7RbkYb8D5JekHRNntebJKmzIuMquyX9VlJjSdl59uQfPVEB64K8Lmm1pP6SjkhqLumM4rwzAJAsaLwBIABmViH6YgNJoyUtlTRLkpxzeyUNz7c3U9JHZnaec2579Iz19ZL6OOdeiu57X9IXkvaXMtuPJDWUdJ1zbnX04XmlOSYAJCJGTQAg/tVS5Kx0tqQNklIl/dI598OZajO7ycyWmtm30X0fR/+pWfTvDoqcuZ517HWccwclzSmDfPskbZM01sz6mNlZZXBMAEg4NN4AEP+yJF2kyJ1Nhko6VdL0Y/9oZr+QNFnSJ5J6Rff9QpJJqhzddo6kb51zR/Ide09pwznnnKRukr6UNEHSLjP70MxSS3tsAEgkjJoAQPw76pxbFn15sZkdljTFzHo7515UpNle4Jy769grmFnnfMfYJamamZ2Sr/k+u4i3fTj69yn5Hq8h6atjC+fcekm9o2Mul0v6L0VuTXhe0e8eACQHzngDQGCcc1MVuZDxvuhDVSR9n2/bQB1/UeRiRc6AX3/sATM7XZEz1YXZHn298/O8Xl1JLU6QLcc5976kRyXVLqN7jQNAQuCMNwCE6a+SpppZV0XmtJ+K3h5woaQekrrm3eycW2NmsyQ9Y2ZnKnIG/DeSDhb2RpxzO8xsiaSHzOyQpAqSfqc89/qOPovmI5JmKHKnlZqK/FCQ4ZzLKot3FgASAY03AIRphqQHFbklYA9F7ioyUpGZ7ncl9ZO0IN/rDJY0RtJjkg5IekqRW/71KuJt9ZU0XtJzipwB/zdJ9+T5913RP7+XVEeRmfS5ku4/qfcMABKURa6JAQAAAFCemPEGAAAAYoDGGwAAAIgBGm8AAAAgBmi8AQAAgBig8QYAAABigMYbAAAAiAEabwAAACAGaLwBAACAGPg/RZ/1FCxR50YAAAAASUVORK5CYII=",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fef36defb10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig = plt.figure(figsize=(12,8))\n",
    "\n",
    "data=np.loadtxt('../postProcessing/sampleDict/20/s2_U.xy', skiprows=0)\n",
    "\n",
    "plt.plot(data[:,0],data[:,1],'o',label='icoFoam')\n",
    "\n",
    "vmax = max(data[:,1])\n",
    "rmax = 0.5\n",
    "\n",
    "x = np.linspace(-1*rmax, rmax, 100)\n",
    "sol = vmax*(1 - x**2/rmax**2)\n",
    "\n",
    "plt.plot(x,sol,'-',label='Analytical solution')\n",
    "\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "plt.xlabel('Radius',fontsize=15)\n",
    "plt.ylabel('$V_{axial}$',fontsize=15)\n",
    "\n",
    "\n",
    "\n",
    "#To save a figure\n",
    "#plt.savefig('test.png', format='png', dpi=300)\n",
    "#plt.savefig('test.pdf', format='pdf', dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For completeness, I will give you a few instructtions of how to work with arrays in numpy.<br>\n",
    "\n",
    "To create an array (in this case a vector):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1, 2, 3, 4])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.array([1, 2, 3, 4])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To create a multi-dimensional array and save the values in the variable x:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x = np.array([[1, 2, 3, 4], [5, 6, 7, 8]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To know the type of a variable you use type():"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "numpy.ndarray"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which states that the variable x is a numpy ndarray. <br>\n",
    "To know the dimensions of the array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2, 4)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To print one value of the array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x[0,0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To print one column of the array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1, 5])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x[:,0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and to print a row:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1, 2, 3, 4])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x[0,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To compute the mean of the previous row:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.5"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(x[0,:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "tttt_playground",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
